\documentclass[12pt]{article}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{fullpage}
\usepackage{amsmath}
\usepackage{amsthm}
\usepackage{array}
\usepackage{color}
\usepackage{graphicx}
\usepackage{float}
\usepackage[hidelinks]{hyperref}
\usepackage{listings}
\usepackage[margin=1in]{geometry}
\usepackage{setspace}
\usepackage{natbib}
\usepackage{proof}
\usepackage{multirow}
\usepackage{hhline}
\usepackage{wrapfig}
\usepackage [english]{babel}
\usepackage{grffile}
\usepackage{tikz}
\usepackage{pgfplots}
\usepackage[tikz]{bclogo}
\usetikzlibrary{chains}
\usetikzlibrary{positioning}
\usetikzlibrary{arrows}
\usepackage{lscape}
\usepackage [autostyle, english = american]{csquotes}
\usepackage{enumitem}
\usepackage{caption}
\captionsetup{justification   = raggedright,
              singlelinecheck = false}
\usepackage{subcaption}
\usepackage{indentfirst}
\usepackage{hyperref}
\usepackage{pdfpages}
\usepackage{booktabs}
\usepackage{pgfgantt}
\newcommand{\tabitem}{~~\llap{\textbullet}~~}
\bibliographystyle{apsr}
\bibpunct{(}{)}{;}{a}{,}{,}
\DeclareGraphicsExtensions{.pdf,.png,.jpg}
\setlength{\tabcolsep}{.18cm}
\usepackage{fancyvrb}
\usepackage{etoc}
\usepackage{sectsty}
\partfont{\fontsize{13}{13}\selectfont}
\sectionfont{\fontsize{13}{13}\selectfont}
\subsectionfont{\fontsize{12}{12}\selectfont}

\usepackage{titlesec}

\makeatletter
\@addtoreset{section}{part}
\@addtoreset{figure}{part}
\@addtoreset{table}{part}
\makeatother
\titleformat{\part}[display]
{\normalfont\LARGE\bfseries}{}{0pt}{}

\makeatletter
\renewcommand{\l@section}{\@dottedtocline{1}{1.5em}{2.6em}}
\renewcommand{\l@subsection}{\@dottedtocline{2}{4.0em}{3.6em}}
\renewcommand{\l@subsubsection}{\@dottedtocline{3}{7.4em}{4.5em}}
\makeatother

% === new commands ===
\newcommand\ud{\mathrm{d}}
\newcommand\dist{\buildrel\rm d\over\sim}
\newcommand\ind{\stackrel{\rm indep.}{\sim}}
\newcommand\iid{\stackrel{\rm i.i.d.}{\sim}}
\newcommand\logit{{\rm logit}}
\renewcommand\r{\right}
\renewcommand\l{\left}
\newcommand\Var{{\rm Var}}
\newcommand\var{{\rm var}}
\newcommand\Cov{{\rm Cov}}
\newcommand\E{\mathbb{E}}
\newcommand\V{\mathbb{V}}
\newcommand\cN{\mathcal{N}}
\newcommand\cS{\mathcal{S}}
\newcommand\cX{\mathcal{X}}
\newcommand\wX{\widetilde{X}}
\newcommand{\diag}{\mathop{\mathrm{diag}}}
\newcommand{\argmax}{\operatornamewithlimits{argmax}}
\newcommand{\argmin}{\operatornamewithlimits{argmin}}
\newcommand{\indep}{\mbox{$\perp\!\!\!\perp$}}
\def\independenT#1#2{\mathrel{\rlap{$#1#2$}\mkern2mu{#1#2}}}
\DeclareMathOperator{\sgn}{sgn}
\providecommand{\norm}[1]{\lVert#1\rVert}

\def\changemargin#1#2{\list{}{\rightmargin#2\leftmargin#1}\item[]}
\let\endchangemargin=\endlist 

\newcolumntype{L}[1]{>{\raggedright\let\newline\\\arraybackslash\hspace{0pt}}m{#1}}
\newcolumntype{C}[1]{>{\centering\let\newline\\\arraybackslash\hspace{0pt}}m{#1}}
\newcolumntype{R}[1]{>{\raggedleft\let\newline\\\arraybackslash\hspace{0pt}}m{#1}}
\newcolumntype{P}[1]{>{\raggedright\tabularxbackslash}p{#1}}

\newcounter{myWeekNum} %for gantt
\stepcounter{myWeekNum}
%
\newcommand{\myWeek}{\themyWeekNum
    \stepcounter{myWeekNum}
    \ifnum\themyWeekNum=53
         \setcounter{myWeekNum}{1}
    \else\fi
}

% ==Cross Referencing Different Docs
\usepackage{xr}
\externaldocument{Colombia_Paper}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title{Left Out: How Political Ideology Affects Support for\\ Migrants in Colombia:\\
Supplementary Information\footnote{This research received institutional review board (IRB) approval from 
Princeton University ($\#$8335), UCLA (\#19-001733), and UBC ($\#$H19-03288). 
Our Pre-Analysis Plan was archived in the OSF repository \url{https://osf.io/v2db9/}. 
All replication material, including {\tt R} code and data, will be made available via Harvard University's Dataverse.
}
}

\author{Alisha Holland\thanks{Associate Professor, Government Department, Harvard University, \href{mailto:aholland@fas.harvard.edu}{aholland@fas.harvard.edu}, \href{http://alishaholland.com/about/}{www.alishaholland.com}}
\hspace{1.5cm}
Margaret E. Peters\thanks{Professor, Department of Political Science, UCLA, \href{mailto:mepeters@ucla.edu}{mepeters@ucla.edu}, \href{http://www.maggiepeters.com/}{www.maggiepeters.com}}
\hspace{1.5cm}
Yang-Yang Zhou\thanks{Assistant Professor, Department of Government, Dartmouth College, \href{mailto:yang-yang.zhou@dartmouth.edu}{yang-yang.zhou@dartmouth.edu}, \href{https://www.yangyangzhou.com/}{www.yangyangzhou.com}}
}

\date{\today}

%%%%%%%%%%%%%%%%% END OF PREAMBLE %%%%%%%%%%%%%%%%

\begin{document}

<<eval=TRUE,echo=FALSE, results='hide', message=FALSE>>= 
library(knitr)

opts_chunk$set(cache = TRUE, 
               cache.path = 'cache_SI/',
               fig.path = 'figures_SI/', 
               tidy = TRUE, 
               echo = FALSE, 
               warning = FALSE, 
               message = FALSE, 
               fig.pos = 'H',
               dev = 'pdf', 
               dpi=200)

options(width = 110, digits = 2)

@


\maketitle

<<eval=TRUE, echo = FALSE, tidy=TRUE, warning=FALSE, error=FALSE, message=FALSE>>=

setwd("Paper_Inputs")

## Load packages and functions
source("functions.R")

## Load cleaned data
df_col <- readRDS("colombia_clean.RDS") 
df_ven <- readRDS("venezuelans_clean.RDS") 

##Additional Cleaning

## Colombians
df_col <- df_col %>% # rescale these variables to min 0 max 1 for observational regression analysis
  mutate(dir_contact_bi = case_when(dir_contact_index > 1~1, 
                                    dir_contact_index <= 1~0,
                                           TRUE~NA_real_),
         ven_beg_bi = case_when(ven_beg > 2~1, 
                                ven_beg <= 2~0,
                                TRUE~NA_real_),
         open_index_res = scales::rescale(open_index), 
         partisanship_res = scales::rescale(partisanship),
         skilled_labor_res = scales::rescale(skilled_labor),
         contract_res = scales::rescale(contract),
         salary_res = scales::rescale(salary),
         benefits_index_res = scales::rescale(benefits_index),
         dir_contact_index_res = scales::rescale(dir_contact_index),
         indir_contact_index_res = scales::rescale(indir_contact_index),
         cultural_index_res = scales::rescale(cultural_index)
         ) 

df_col_cali <- df_col[df_col$city == "Cali",]
df_col_cucuta <- df_col[df_col$city == "Cúcuta",]

df_col_left <- df_col[df_col$partisanship == 1,]
df_col_right <- df_col[df_col$partisanship != 1,]

df_col_contact <- df_col[df_col$dir_contact_bi == 1,]
df_col_nocontact <- df_col[df_col$dir_contact_bi != 1,]

## Venezuelans
df_ven<-df_ven %>%
  mutate(left_partisanship = case_when(partisanship == 1~1,
                                partisanship >1~0,
                                TRUE~NA_real_)
         )

df_ven_cali <- df_ven[df_ven$city == "Cali",]
df_ven_cucuta <- df_ven[df_ven$city == "Cúcuta",]

## Load LAPOP merged data
col.merge <- readRDS("col_merged.RDS") # Colombia
ven.merge <- readRDS("ven_merged.RDS") # Venezuela

@

%\newpage
%To number supplemental material with 'S': 
\renewcommand{\thesection}{S\arabic{section}}   
\renewcommand{\thetable}{S\arabic{table}}   
\renewcommand{\thefigure}{S\arabic{figure}}
\renewcommand{\theequation}{S\arabic{equation}}

%\part*{\large For Online Publication: Supplementary Information}

\addtocontents{toc}{\protect\setcounter{tocdepth}{2}}
\tableofcontents

\pagenumbering{gobble}

\newpage
\pagenumbering{arabic}
\setcounter{page}{1}
\setstretch{1}

\newpage
\section{Examples of the Importance of Migrants' Political Views in Multiple Contexts}
\label{SIsec:examples}

There are many historical examples in which political concerns have affected the treatment of immigrants. The U.S., for instance, has a long history of limiting immigration and citizenship rights in reaction to political fears. In the colonial and early Republic period of the U.S., concerns that Catholic immigrants would be loyal to the Pope, rather than to the English monarch and later to the Republic, led to laws prohibited Catholics from naturalization or serving as elected officials under the New York Constitution unless they renounced their faith.  These laws stayed in place until 1806 \citep{duncan2005citizens}. In the wake of the French Revolution, the arrival of refugees from France and radical sympathizers from Great Britain and Ireland led to fears that immigrants would spread radical ideas, leading to the passage of the Alien and Sedition Acts \citep[p.662]{cogliano1999america}. 

Political concerns about immigrants resurfaced in the twentieth century around the spread of socialist, communist, and anarchist ideas \citep{higham1983strangers}. In Europe, a spate of political assassinations by anarchists led to deportations and increased political vetting of immigrants. In the U.S., the Alien Exclusion Act of 1903 banned anarchists and a 1906 law denaturalized anarchists \citep[p.59]{kraut2020threat}. At the height of the Red Scare of the 1920s, approximately 3,000 immigrants were held as radicals at Ellis Island and 556 were deported \citep[p.74]{kraut2020threat} and concerns about the spread of leftist ideas help motivate the 1921 and 1924 Quota Acts \citep{kraut2012global, kraut2020threat}. 

Political fears can also involve misperceptions about the views of entire national groups. A clear example  comes from the internment of Japanese-Americans as potential fascists, even though most had come decades before or were born in the U.S. Similarly, Vietnamese refugees were often portrayed as communist enemies in the 1970s, even though most were fleeing the policies of their government \citep[e.g.][]{wooten1975}. 

In the developing world, migrants often leave countries pursuing extreme ideological projects and cross to neighboring countries. Host countries may worry about possible ``contagion'' in which political movements in a neighboring country spread to their own. For instance, in Thailand the government and media portrayed Vietnamese refugees as communist sympathizers and a potential ``vanguard'' for a Communist invasion, even while the majority of refugees disavowed their government, as a way to discredit local socialist groups \citep[39]{flood1977vietnamese}. 

Work on emigration begins from the opposite assumption, that migrants are critics of their home governments (for a review, see \citet{kapur2014political}). Governments often consider emigration a useful safety valve to relieve political tensions and the exit of political opponents--and sometimes the deliberate exile of high-profile dissidents--can stabilize authoritarian governments \citep{miller2020restraining}. From this perspective, migrants are more likely to be political dissidents or at least come to oppose the policies that force them to leave. 


\newpage
\section{Respondent Demographic Covariates}
\label{SIsec:demog}

\subsection{Respondent Summary Statistics}

<<Sumstats_Col, eval = TRUE, echo=FALSE, fig.align='center', warning=FALSE, message=FALSE, strip.white=TRUE, results='asis'>>=

descriptlabels <-  c("Male",
                       "Age Group",
                       "Race", 
                       "Education", 
                       "Number of Children", 
                       "Marriage Status",
                       "Religion", 
                       "Religiosity",
                       "Wealth Index",
                       "Labor Contract", 
                       "Salary",
                       "Political Ideology",
                       "Public Benefits Index",
                       "Indirect Contact Index",
                       "Direct Contact Index",
                       "Internally Displaced",
                       "Family Internally Displaced")

col_descript <- subset(df_col, 
            select = c("city", 
                       "male",
                       "age",
                       "race", 
                       "education", 
                       "kids", 
                       "marriage",
                       "religion2", 
                       "religiosity",
                       "wealth_index",
                       "contract",
                       "salary",
                       "partisanship",
                       "benefits_index",
                       "dir_contact_index",
                       "indir_contact_index",
                       "int_disp", 
                       "int_disp_fam"))

stargazer(col_descript,
          title = "Descriptive statistics of Colombian respondents",
          median = TRUE,
          digits=2, 
          digits.extra=2,
          covariate.labels = descriptlabels,
          font.size = "small",
          label = "tab:sumstats_col")

@

<<Sumstats_Ven, eval = TRUE, echo=FALSE, fig.align='center', warning=FALSE, message=FALSE, strip.white=TRUE, results='asis'>>=

descriptlabels <-  c("Male",
                       "Age Group",
                       "Race", 
                       "Education", 
                       "Number of Children", 
                       "Marriage Status",
                       "Religion", 
                       "Religiosity",
                       "Wealth Index",
                       "Labor Contract", 
                       "Salary",
                       "Political Ideology",
                       "Public Benefits Index")

col_descript <- subset(df_ven, 
            select = c("city", 
                       "male",
                       "age",
                       "race", 
                       "education", 
                       "kids", 
                       "marriage",
                       "religion2", 
                       "religiosity",
                       "wealth_index_col",
                       "contract",
                       "salary",
                       "partisanship",
                       "benefits_index"))

stargazer(col_descript,
          title = "Descriptive statistics of Venezuelan respondents",
          median = TRUE,
          digits=2, 
          digits.extra=2,
          covariate.labels = descriptlabels,
          font.size = "small",
          label = "tab:sumstats_ven")

@


\newpage
\subsection{Demographics for Colombians}
<<demographics_colombians, eval = TRUE, echo = FALSE, tidy=TRUE, fig.width = 16, fig.height = 18, out.width= "1\\linewidth", fig.align='center', warning=FALSE, message=FALSE, strip.white=TRUE, fig.cap="Colombian respondent demographics.">>=


# City
city_col <- ggplot(df_col[is.na(df_col$city)==FALSE,],
                   aes(x=city)) +
        geom_bar(aes(y = (..count..)/sum(..count..))) +
        ggtitle("City") +
        #ylim(0, .6) +
        xlab("") +
        ylab("Proportion") +
        #scale_fill_manual(values=c("white", "gray60")) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Birth Country (all Colombia)
#table(df_col$birth_country)

 # Gender
gender_col <- ggplot(df_col[is.na(df_col$male)==FALSE,],
                   aes(x=male)) +
        geom_bar(aes(y = (..count..)/sum(..count..))) +
        ggtitle("Gender") +
        #ylim(0,1) +
        xlab("") +
        ylab("Proportion") +
        #scale_fill_manual(values=c("white", "gray60")) +
        scale_x_continuous(breaks=seq(0, 1, 1), 
                  labels=c("Female", 
                           "Male")) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Age
age_col <- ggplot(df_col[is.na(df_col$age)==FALSE,],
                   aes(x=age)) +
        geom_bar(aes(y = (..count..)/sum(..count..))) +
        ggtitle("Age") +
        #ylim(0,1) +
        xlab("") +
        ylab("Proportion") +
        #scale_fill_manual(values=c("white", "gray60")) +
        scale_x_continuous(breaks=seq(1, 5, 1), 
                  labels=c("18-24", 
                           "25-34",
                           "35-44",
                           "45-54",
                           "55+")) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Race
race_col <- ggplot(df_col[is.na(df_col$race)==FALSE,],
                   aes(x=race)) +
        geom_bar(aes(y = (..count..)/sum(..count..))) +
        ggtitle("Race") +
        #ylim(0,1) +
        xlab("") +
        ylab("Proportion") +
        #scale_fill_manual(values=c("white", "gray60")) +
        scale_x_continuous(breaks=seq(1, 6, 1), 
                  labels=c("White", 
                           "Mestizo",
                           "Indigenous",
                           "Afro/\nBlack",
                           "Mulato",
                           "Other")) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Education
edu_col <- ggplot(df_col[is.na(df_col$education)==FALSE,],
                   aes(x=education)) +
        geom_bar(aes(y = (..count..)/sum(..count..))) +
        ggtitle("Education") +
        #ylim(0,1) +
        xlab("") +
        ylab("Proportion") +
        #scale_fill_manual(values=c("white", "gray60")) +
        scale_x_continuous(breaks=seq(0, 5, 1), 
                  labels=c("None",
                           "Pri", 
                           "Sec",
                           "Tech",
                           "Univ",
                           "Masters/\nPhD")) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Number of Children
kids_col <- ggplot(df_col, aes(x=kids)) +
        geom_histogram(aes(y=..density..), binwidth=1, colour="white") +
        #geom_density(alpha=.2) +
        ggtitle("Number of Children") +
        xlab("") +
        ylab("Density") +
        #ylim(0, .15) +
        #scale_x_continuous(breaks=seq(10, 50, 10)) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Marriage
marriage_col <- ggplot(df_col[is.na(df_col$marriage)==FALSE,],
                   aes(x=marriage)) +
        geom_bar(aes(y = (..count..)/sum(..count..))) +
        ggtitle("Marriage Status") +
        #ylim(0,1) +
        xlab("") +
        ylab("Proportion") +
        #scale_fill_manual(values=c("white", "gray60")) +
        scale_x_continuous(breaks=seq(1, 5, 1), 
                  labels=c("Single",
                           "Married", 
                           "Divorced",
                           "Civil\nUnion",
                           "Widow")) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Religion
religion_col <- ggplot(df_col[is.na(df_col$religion2)==FALSE,],
                   aes(x=religion2)) +
        geom_bar(aes(y = (..count..)/sum(..count..))) +
        ggtitle("Religion") +
        #ylim(0,1) +
        xlab("") +
        ylab("Proportion") +
        #scale_fill_manual(values=c("white", "gray60")) +
        scale_x_continuous(breaks=seq(1, 3, 1), 
                  labels=c("Catholic",
                           "Evangelical", 
                           "Other")) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Religiosity
religiosity_col <- ggplot(df_col[is.na(df_col$religiosity)==FALSE,],
                   aes(x=religiosity)) +
        geom_bar(aes(y = (..count..)/sum(..count..))) +
        ggtitle("Religiosity") +
        #ylim(0,1) +
        xlab("") +
        ylab("Proportion") +
        #scale_fill_manual(values=c("white", "gray60")) +
        scale_x_continuous(breaks=seq(0, 3, 1), 
                  labels=c("Unimportant",
                           "Somewhat\nunimportant", 
                           "Somewhat\nimportant",
                           "Very\nimportant")) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Wealth Index
wealth_col <- ggplot(df_col, aes(x=wealth_index)) +
        geom_histogram(aes(y=..density..), binwidth=1, colour="white") +
        #geom_density(alpha=.2) +
        ggtitle("Wealth Index") +
        xlab("") +
        ylab("Density") +
        #ylim(0, .15) +
        scale_x_continuous(breaks=seq(0, 13, 5)) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Occupation skill
skill_col <- ggplot(df_col[is.na(df_col$skilled_labor)==FALSE,],
                   aes(x=skilled_labor)) +
        geom_bar(aes(y = (..count..)/sum(..count..))) +
        ggtitle("Labor skill") +
        #ylim(0,1) +
        xlab("") +
        ylab("Proportion") +
        #scale_fill_manual(values=c("white", "gray60")) +
        scale_x_continuous(breaks=seq(0, 2, 1), 
                  labels=c("Unemployed",
                           "Low skilled", 
                           "High skilled")) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Contract
contract_col <- ggplot(df_col[is.na(df_col$contract)==FALSE,],
                   aes(x=contract)) +
        geom_bar(aes(y = (..count..)/sum(..count..))) +
        ggtitle("Labor Contract") +
        #ylim(0,1) +
        xlab("") +
        ylab("Proportion") +
        #scale_fill_manual(values=c("white", "gray60")) +
        scale_x_continuous(breaks=seq(1, 3, 1), 
                  labels=c("Informal",
                           "Formal, fixed", 
                           "Formal, temp")) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Salary
salary_col <- ggplot(df_col[is.na(df_col$salary)==FALSE,],
                   aes(x=salary)) +
        geom_bar(aes(y = (..count..)/sum(..count..))) +
        ggtitle("Salary") +
        #ylim(0,1) +
        xlab("") +
        ylab("Proportion") +
        #scale_fill_manual(values=c("white", "gray60")) +
        scale_x_continuous(breaks=seq(1, 4, 1), 
                  labels=c("Not enough,\nMajor\ndifficulties",
                           "Not enough,\nSome\ndifficulties",
                           "Enough,\nSome\ndifficulties",
                           "Enough,\nSavings")) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Partisanship
partisanship_col <- ggplot(df_col[is.na(df_col$partisanship)==FALSE,],
                   aes(x=partisanship)) +
        geom_bar(aes(y = (..count..)/sum(..count..))) +
        ggtitle("Partisanship") +
        #ylim(0,1) +
        xlab("") +
        ylab("Proportion") +
        #scale_fill_manual(values=c("white", "gray60")) +
        scale_x_continuous(breaks=seq(1, 3, 1), 
                  labels=c("Left",
                           "Center",
                           "Right")) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Presidential Vote
## something wrong with responses, on scale from 1 to 5

# Benefits
benefits_col <- ggplot(df_col, aes(x=benefits_index)) +
        geom_histogram(aes(y=..density..), binwidth=1, colour="white") +
        #geom_density(alpha=.2) +
        ggtitle("Number of Public Benefits") +
        xlab("") +
        ylab("Density") +
        #ylim(0, .15) +
        scale_x_continuous(breaks=seq(0, 5, 1)) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Direct Contact 
directcontact_col <- ggplot(df_col, aes(x=dir_contact_index)) +
        geom_histogram(aes(y=..density..), binwidth=.5, colour="white") +
        #geom_density(alpha=.2) +
        ggtitle("Number of Venezuelan Friends\nor Colleagues (Direct Contact)") +
        xlab("") +
        ylab("Density") +
        #ylim(0, .15) +
        scale_x_continuous(breaks=seq(0, 3, 1), 
                  labels=c("0",
                           "2",
                           "4",
                           "6+")) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Indirect contact
indirectcontact_col <- ggplot(df_col, aes(x=indir_contact_index)) +
        geom_histogram(aes(y=..density..), binwidth=.5, colour="white") +
        #geom_density(alpha=.2) +
        ggtitle("Frequency of Seeing\nVenezuelans (Indirect Contact)") +
        xlab("") +
        ylab("Density") +
        #ylim(0, .15) +
        scale_x_continuous(breaks=seq(0, 3, 1), 
                  labels=c("Never",
                           "Sometimes",
                           "Frequently",
                           "All the\ntime")) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )


## Plot
city_col + gender_col + age_col + race_col + 
  edu_col + kids_col + marriage_col + religion_col + 
  religiosity_col + wealth_col + contract_col + salary_col +
  partisanship_col + benefits_col + directcontact_col + indirectcontact_col +
  plot_layout(ncol=4)


@

\newpage
\subsection{Demographics for Venezuelans}
<<demographics_venezuelans, eval = TRUE, echo = FALSE, tidy=TRUE, fig.width = 16, fig.height = 18, out.width= "1\\linewidth", fig.align='center', warning=FALSE, message=FALSE, strip.white=TRUE, fig.cap="Venezuelan respondent demographics.">>=

# City
city_ven <- ggplot(df_ven[is.na(df_ven$city)==FALSE,],
                   aes(x=city)) +
        geom_bar(aes(y = (..count..)/sum(..count..))) +
        ggtitle("City") +
        #ylim(0, .6) +
        xlab("") +
        ylab("Proportion") +
        #scale_fill_manual(values=c("white", "gray60")) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Birth Country (all Venezuela)
#table(df_ven$birth_country)

 # Gender
gender_ven <- ggplot(df_ven[is.na(df_ven$male)==FALSE,],
                   aes(x=male)) +
        geom_bar(aes(y = (..count..)/sum(..count..))) +
        ggtitle("Gender") +
        #ylim(0,1) +
        xlab("") +
        ylab("Proportion") +
        #scale_fill_manual(values=c("white", "gray60")) +
        scale_x_continuous(breaks=seq(0, 1, 1), 
                  labels=c("Female", 
                           "Male")) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Age
age_ven <- ggplot(df_ven[is.na(df_ven$age)==FALSE,],
                   aes(x=age)) +
        geom_bar(aes(y = (..count..)/sum(..count..))) +
        ggtitle("Age") +
        #ylim(0,1) +
        xlab("") +
        ylab("Proportion") +
        #scale_fill_manual(values=c("white", "gray60")) +
        scale_x_continuous(breaks=seq(1, 5, 1), 
                  labels=c("18-24", 
                           "25-34",
                           "35-44",
                           "45-54",
                           "55+")) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Race
race_ven <- ggplot(df_ven[is.na(df_ven$race)==FALSE,],
                   aes(x=race)) +
        geom_bar(aes(y = (..count..)/sum(..count..))) +
        ggtitle("Race") +
        #ylim(0,1) +
        xlab("") +
        ylab("Proportion") +
        #scale_fill_manual(values=c("white", "gray60")) +
        scale_x_continuous(breaks=seq(1, 6, 1), 
                  labels=c("White", 
                           "Mestizo",
                           "Indigenous",
                           "Afro/\nBlack",
                           "Mulato",
                           "Other")) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Education
edu_ven <- ggplot(df_ven[is.na(df_ven$education)==FALSE,],
                   aes(x=education)) +
        geom_bar(aes(y = (..count..)/sum(..count..))) +
        ggtitle("Education") +
        #ylim(0,1) +
        xlab("") +
        ylab("Proportion") +
        #scale_fill_manual(values=c("white", "gray60")) +
        scale_x_continuous(breaks=seq(0, 5, 1), 
                  labels=c("None",
                           "Pri", 
                           "Sec",
                           "Tech",
                           "Univ",
                           "Masters/\nPhD")) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Number of Children
kids_ven <- ggplot(df_ven, aes(x=kids)) +
        geom_histogram(aes(y=..density..), binwidth=1, colour="white") +
        #geom_density(alpha=.2) +
        ggtitle("Number of Children") +
        xlab("") +
        ylab("Density") +
        #ylim(0, .15) +
        #scale_x_continuous(breaks=seq(10, 50, 10)) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Marriage
marriage_ven <- ggplot(df_ven[is.na(df_ven$marriage)==FALSE,],
                   aes(x=marriage)) +
        geom_bar(aes(y = (..count..)/sum(..count..))) +
        ggtitle("Marriage Status") +
        #ylim(0,1) +
        xlab("") +
        ylab("Proportion") +
        #scale_fill_manual(values=c("white", "gray60")) +
        scale_x_continuous(breaks=seq(1, 5, 1), 
                  labels=c("Single",
                           "Married", 
                           "Divorced",
                           "Civil\nUnion",
                           "Widow")) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Religion
religion_ven <- ggplot(df_ven[is.na(df_ven$religion2)==FALSE,],
                   aes(x=religion2)) +
        geom_bar(aes(y = (..count..)/sum(..count..))) +
        ggtitle("Religion") +
        #ylim(0,1) +
        xlab("") +
        ylab("Proportion") +
        #scale_fill_manual(values=c("white", "gray60")) +
        scale_x_continuous(breaks=seq(1, 3, 1), 
                  labels=c("Catholic",
                           "Evangelical", 
                           "Other")) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Religiosity
religiosity_ven <- ggplot(df_ven[is.na(df_ven$religiosity)==FALSE,],
                   aes(x=religiosity)) +
        geom_bar(aes(y = (..count..)/sum(..count..))) +
        ggtitle("Religiosity") +
        #ylim(0,1) +
        xlab("") +
        ylab("Proportion") +
        #scale_fill_manual(values=c("white", "gray60")) +
        scale_x_continuous(breaks=seq(0, 3, 1), 
                  labels=c("Unimportant",
                           "Somewhat\nunimportant", 
                           "Somewhat\nimportant",
                           "Very\nimportant")) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Wealth Index
wealth_ven <- ggplot(df_ven, aes(x=wealth_index_col)) +
        geom_histogram(aes(y=..density..), binwidth=1, colour="white") +
        #geom_density(alpha=.2) +
        ggtitle("Wealth Index (in Colombia)") +
        xlab("") +
        ylab("Density") +
        #ylim(0, .15) +
        scale_x_continuous(breaks=seq(0, 6, 1)) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Occupation skill
skill_ven <- ggplot(df_ven[is.na(df_ven$skilled_labor)==FALSE,],
                   aes(x=skilled_labor)) +
        geom_bar(aes(y = (..count..)/sum(..count..))) +
        ggtitle("Labor skill") +
        #ylim(0,1) +
        xlab("") +
        ylab("Proportion") +
        #scale_fill_manual(values=c("white", "gray60")) +
        scale_x_continuous(breaks=seq(0, 2, 1), 
                  labels=c("Unemployed",
                           "Low skilled", 
                           "High skilled")) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Contract
contract_ven <- ggplot(df_ven[is.na(df_ven$contract)==FALSE,],
                   aes(x=contract)) +
        geom_bar(aes(y = (..count..)/sum(..count..))) +
        ggtitle("Labor Contract") +
        #ylim(0,1) +
        xlab("") +
        ylab("Proportion") +
        #scale_fill_manual(values=c("white", "gray60")) +
        scale_x_continuous(breaks=seq(1, 3, 1), 
                  labels=c("Informal",
                           "Formal, fixed", 
                           "Formal, temp")) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Salary
salary_ven <- ggplot(df_ven[is.na(df_ven$salary)==FALSE,],
                   aes(x=salary)) +
        geom_bar(aes(y = (..count..)/sum(..count..))) +
        ggtitle("Salary") +
        #ylim(0,1) +
        xlab("") +
        ylab("Proportion") +
        #scale_fill_manual(values=c("white", "gray60")) +
        scale_x_continuous(breaks=seq(1, 4, 1), 
                  labels=c("Not enough,\nMajor\ndifficulties",
                           "Not enough,\nSome\ndifficulties",
                           "Enough,\nSome\ndifficulties",
                           "Enough,\nSavings")) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Partisanship
partisanship_ven <- ggplot(df_ven[is.na(df_ven$partisanship)==FALSE,],
                   aes(x=partisanship)) +
        geom_bar(aes(y = (..count..)/sum(..count..))) +
        ggtitle("Partisanship") +
        #ylim(0,1) +
        xlab("") +
        ylab("Proportion") +
        #scale_fill_manual(values=c("white", "gray60")) +
        scale_x_continuous(breaks=seq(1, 3, 1), 
                  labels=c("Left",
                           "Center",
                           "Right")) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Presidential Vote
## something wrong with responses, on scale from 1 to 5

# Benefits
benefits_ven <- ggplot(df_ven, aes(x=benefits_index)) +
        geom_histogram(aes(y=..density..), binwidth=1, colour="white") +
        #geom_density(alpha=.2) +
        ggtitle("Number of Public Benefits") +
        xlab("") +
        ylab("Density") +
        #ylim(0, .15) +
        scale_x_continuous(breaks=seq(0, 5, 1)) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

# Main Reason for Leaving Venezuela
reasonleave_ven <- ggplot(df_ven[is.na(df_ven$reason_leave_num1)==FALSE & 
                                   df_ven$reason_leave_num1 < 11,],
                   aes(x=reason_leave_num1)) +
        geom_bar(aes(y = (..count..)/sum(..count..))) +
        ggtitle("Main Reason for Leaving Venezuela") +
        #ylim(0,1) +
        xlab("") +
        ylab("Proportion") +
        #scale_fill_manual(values=c("white", "gray60")) +
        scale_x_continuous(breaks=seq(1, 10, 1), 
                  labels=c("Economic",
                           "Security",
                           "Unemployment",
                           "Food",
                           "Medicine",
                           "Better life",
                           "Join family",
                           "Pregnant",
                           "Help family",
                           "Education")) +
        theme(
        panel.background = element_blank(),
        legend.position = "none",
        panel.border = element_rect(colour = "gray", fill=NA, size=.8)
        )

## Plot
(city_ven + gender_ven + age_ven + race_ven + plot_layout(ncol=4)) / 
(edu_ven + kids_ven + marriage_ven + religion_ven + plot_layout(ncol=4)) /
(religiosity_ven + wealth_ven + contract_ven + salary_ven + plot_layout(ncol=4)) /
(partisanship_ven + benefits_ven + reasonleave_ven + plot_layout(ncol=3,widths=c(1,1,2))) 

@


\newpage
\section{Sample Representativeness: Comparing Respondents with LAPOP}
\label{SIsec:lapop}

\subsection{Comparison for Colombians}

We use LAPOP Colombia as a comparison, which is a self-weighted, nationally representative sample with 1,663 respondents. For the salary comparison, we use LAPOP (2016), and all other variables are from LAPOP (2018). Compared to a national sample of Colombians, our Colombian respondents tended to have less children, be slightly more politically left (although the average is still towards the right), and wealthier -- reflecting that our respondents are urban. They were also slightly older, less likely to be a student or employed, and more likely to be Christian (Catholic or Protestant).

<<lapopbalance_colombians, eval = TRUE, echo = FALSE, tidy=TRUE>>=

## Make Colombia Table 

###set vars as either numeric or factor
col.merge <- col.merge %>% 
  mutate_at(c("male", "kids", "pol_ideo", 
              "wealth_car", "wealth_washing", "wealth_tv",  
              "wealth_computer", "wealth_cell",  
              "wealth_internet","contract"),~as.numeric(as.character(.x)))

col.merge <- col.merge %>% 
  mutate_at(c("age", "marriage", "education", "race", 
              "salary", "religion2", "religiosity", "job"),~as.factor(.x))
    
## Create long dataset - numeric
col.merge.num <- col.merge %>% 
  dplyr::select(male, kids, pol_ideo, wealth_car, wealth_washing, 
                wealth_tv, wealth_computer, wealth_cell, 
                wealth_internet, contract, lapop) %>%
  pivot_longer(cols = -lapop, names_to = "variable", values_to = "value") %>%
  drop_na()

bal.tab.num.col <- col.merge.num %>%
  group_by(variable) %>%
  do({
    
    ttest_out <- t.test(value ~ lapop, data = .)
    n_lapop <- sum(.$lapop)
    n_hpz <- nrow(.) - n_lapop
    data.frame(mean_lapop = ttest_out$estimate[1], 
               #se_lapop = sqrt(var(.$value[.$lapop == 1]) / n_lapop),
               mean_hpz = ttest_out$estimate[2],
               #se_hpz = sqrt(var(.$value[.$lapop == 0]) / n_hpz),
               diff_means = ttest_out$estimate[1] - ttest_out$estimate[2],
               pval_difference = ttest_out$p.value)
}) %>%
  rename(`Variable` = variable, 
         `LAPOP Mean` = mean_lapop, 
         #`LAPOP SE` = se_lapop, 
         `Sample Mean` = mean_hpz,
         #`Sample SE` = se_hpz,
         `Diff. in Means` = diff_means,
         `P-Value` = pval_difference) %>%
  mutate(Variable = case_when(Variable == "contract"~"Formal Contract",
                              Variable == "kids"~"Number of Children",
                              Variable == "male"~"Male",
                              Variable == "pol_ideo"~"Political Ideology (1 left - 10 right)",
                              Variable == "wealth_car"~"Owns Car",
                              Variable == "wealth_cell"~"Owns Cellphone",
                              Variable == "wealth_computer"~"Owns Computer",
                              Variable == "wealth_internet"~"Access to Internet",
                              Variable == "wealth_tv"~"Owns TV",
                              Variable == "wealth_washing"~"Owns Washing Machine"
                              )) %>%
  kable(align = "lccc",
        caption = "Comparing LAPOP and our Sample in Colombia: Numeric Variables",
        label ="tab:lapop_col_num") %>%
  kable_styling(bootstrap_options = "striped")


## Add Variable Labels
col.merge <-col.merge %>% mutate(age = factor(age,
                                             levels = 1:5, 
                                             labels = c("18-24", 
                                                        "25-34",  
                                                        "35-44", 
                                                        "45-54", 
                                                        "55 or above")),
                                marriage   = factor(marriage  , 
                                                    levels = 1:5, 
                                                    labels = c("Single",  
                                                               "Married",
                                                               "Separated/Divorced",
                                                               "Widow",
                                                               "Civil Union")),
                                education   = factor(education, 
                                                     levels = 0:4, 
                                                     labels = c("None",
                                                                "Primary",
                                                                "High School", 
                                                                "University/Technical", 
                                                                "Masters/PhD")),
                                race   = factor(race, 
                                                levels = 1:6, 
                                                labels = c("White",
                                                           "Mestizo",
                                                           "Indigenous", 
                                                           "Black", 
                                                           "Mulatto",
                                                           "Other")),
                                salary   = factor(salary, 
                                                  levels = 1:4, 
                                                  labels = c("Not Enough, Major Difficulties",
                                                             "Not Enough, Some Difficulties" ,
                                                             "Enough Without Great Difficulty",
                                                             "Enough and Can Save")),
                                religion2  = factor(religion2, 
                                                    levels = 1:3, 
                                                    labels = c("Catholic",
                                                               "Protestant" ,
                                                               "Other")),
                                religiosity  = factor(religiosity, 
                                                      levels = 0:3, 
                                                      labels = c("Unimportant",
                                                                 "Somewhat Unimportant",
                                                                 "Somewhat Important",
                                                                 "Very Important")),
                                job  = factor(job, 
                                              levels = 1:4, 
                                              labels = c("Not Seeking Work",
                                                         "Looking for Work",
                                                         "Student",
                                                         "Employed"))
                                
)


## Create long dataset - categorical
col.merge.cat <- col.merge %>% 
  dplyr::select(lapop, age, marriage, education, race, salary, religion2, religiosity, job) %>%
  pivot_longer(cols = -lapop, names_to = "variable", values_to = "value") %>%
  drop_na()

bal.tab.cat.col <- col.merge.cat %>%
  group_by(variable) %>%
  do({
    
    unique_vals <- unique(.$value)
    pval_out <- map_df(unique_vals, function(x){
      x_lapop <- sum(.$value[.$lapop == 1] == x)
      x_hpz <- sum(.$value[.$lapop == 0] == x)
      n_lapop <- sum(.$lapop)
      n_hpz <- nrow(.) - n_lapop
      ptest_out <- prop.test(x = c(x_lapop, x_hpz), n = c(n_lapop, n_hpz))
      data.frame(value = x, 
                 mean_lapop = ptest_out$estimate[1], 
                 #se_lapop = sqrt(ptest_out$estimate[1] * (1 - ptest_out$estimate[1]) / n_lapop),
                 mean_hpz = ptest_out$estimate[2], 
                 #se_hpz = sqrt(ptest_out$estimate[2] * (1 - ptest_out$estimate[2]) / n_hpz),
                 diff_means = ptest_out$estimate[1] - ptest_out$estimate[2],
                 pval_difference = ptest_out$p.value) %>%
        mutate(value = as.factor(value)) 
    }) %>%
      arrange(value)
    
    pval_out  
})  %>%
  rename(`Variable` = variable, 
         `Category` = value, 
         `LAPOP Mean` = mean_lapop, 
         #`LAPOP SE` = se_lapop, 
         `Sample Mean` = mean_hpz,
         #`Sample SE` = se_hpz,
         `Diff. in Means` = diff_means,
         `P-Value` = pval_difference) %>%
  mutate(Variable = case_when(Variable == "age"~"Age",
                              Variable == "education"~"Education",
                              Variable == "job"~"Employment",
                              Variable == "marriage"~"Marriage",
                              Variable == "race"~"Race",
                              Variable == "religion2"~"Religion",
                              Variable == "religiosity"~"Religiosity",
                              Variable == "salary"~"Salary"                              
                              )) %>%
  kable(align = "lccc",
        caption = "Comparing LAPOP and our Sample in Colombia: Categorical Variables",
        label ="tab:lapop_col_cat") %>%
  kable_styling(bootstrap_options = "striped")

# show tables
bal.tab.num.col %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 9)

bal.tab.cat.col %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 9)

@


\newpage
\subsection{Comparison for Venezuelans}

We use LAPOP Venezuela (2016-2017) as a comparison, which is a self-weighted, nationally representative sample with 1,558 respondents. Compared to a national sample of Venezuelans in Venezuela, our respondents of Venezuelan migrants in Colombia tended to be slightly more politically left (although the average is still towards the right), wealthier, younger, more likely to be looking for work and having trouble making ends meet, more Mestizo, less Catholic and more Protestant.

<<lapopbalance_venezuelans, eval = TRUE, echo = FALSE, tidy=TRUE>>=

######     Venezuela Table   #####
###set vars as either numeric or factor
ven.merge <- ven.merge %>% 
  mutate_at(c("male", "kids", "pol_ideo", 
              "wealth_car_ven", "wealth_washing_ven", "wealth_tv_ven",  
              "wealth_computer_ven", "wealth_cell_ven",  
              "wealth_internet_ven"),~as.numeric(as.character(.x)))

ven.merge <- ven.merge %>% 
  mutate_at(c("age", "marriage", "education", "race", 
              "salary", "religion2", "religiosity", "job"),~as.factor(.x))

## Create long dataset - numeric
ven.merge.num <- ven.merge %>% 
  dplyr::select(male, kids, pol_ideo, wealth_car_ven, wealth_washing_ven, 
                wealth_tv_ven, wealth_computer_ven, wealth_cell_ven, 
                wealth_internet_ven, lapop) %>%
  pivot_longer(cols = -lapop, names_to = "variable", values_to = "value") %>%
  drop_na()

bal.tab.num.ven <- ven.merge.num %>%
  group_by(variable) %>%
  do({
    
    ttest_out <- t.test(value ~ lapop, data = .)
    n_lapop <- sum(.$lapop)
    n_hpz <- nrow(.) - n_lapop
    data.frame(mean_lapop = ttest_out$estimate[1], 
               #se_lapop = sqrt(var(.$value[.$lapop == 1]) / n_lapop),
               mean_hpz = ttest_out$estimate[2],
               #se_hpz = sqrt(var(.$value[.$lapop == 0]) / n_hpz),
               diff_means = ttest_out$estimate[1] - ttest_out$estimate[2],
               pval_difference = ttest_out$p.value)
  }) %>%
  rename(`Variable` = variable, 
         `LAPOP Mean` = mean_lapop, 
         #`LAPOP SE` = se_lapop, 
         `Sample Mean` = mean_hpz,
         #`Sample SE` = se_hpz,
         `Diff. in Means` = diff_means,
         `P-Value` = pval_difference) %>%
  mutate(Variable = case_when(Variable == "kids"~"Number of Children",
                              Variable == "male"~"Male",
                              Variable == "pol_ideo"~"Political Ideology (1 left - 10 right)",
                              Variable == "wealth_car_ven"~"Owns Car",
                              Variable == "wealth_cell_ven"~"Owns Cellphone",
                              Variable == "wealth_computer_ven"~"Owns Computer",
                              Variable == "wealth_internet_ven"~"Access to Internet",
                              Variable == "wealth_tv_ven"~"Owns TV",
                              Variable == "wealth_washing_ven"~"Owns Washing Machine"
                              )) %>%
  kable(align = "lccc",
        caption = "Comparing LAPOP and our Sample in Venezuela: Numeric Variables",
        label ="tab:lapop_ven_num") %>%
  kable_styling(bootstrap_options = "striped")  


## Add Variable Labels

ven.merge<-ven.merge %>% mutate(age = factor(age,
                                             levels = 1:5, 
                                             labels = c("18-24", 
                                                        "25-34",  
                                                        "35-44", 
                                                        "45-54", 
                                                        "55 or above")),
                                marriage   = factor(marriage  , 
                                                    levels = 1:5, 
                                                    labels = c("Single",  
                                                               "Married",
                                                               "Separated/Divorced",
                                                               "Widow",
                                                               "Civil Union")),
                                education   = factor(education, 
                                                     levels = 0:4, 
                                                     labels = c("None",
                                                                "Primary",
                                                                "High School", 
                                                                "University/Technical", 
                                                                "Masters/PhD")),
                                race   = factor(race, 
                                                levels = 1:6, 
                                                labels = c("White",
                                                           "Mestizo",
                                                           "Indigenous", 
                                                           "Black", 
                                                           "Mulatto",
                                                           "Other")),
                                salary   = factor(salary, 
                                                  levels = 1:4, 
                                                  labels = c("Not Enough, Major Difficulties",
                                                             "Not Enough, Some Difficulties" ,
                                                             "Enough Without Great Difficulty",
                                                             "Enough and Can Save")),
                                religion2  = factor(religion2, 
                                                    levels = 1:3, 
                                                    labels = c("Catholic",
                                                               "Protestant" ,
                                                               "Other")),
                                religiosity  = factor(religiosity, 
                                                      levels = 0:3, 
                                                      labels = c("Unimportant",
                                                                 "Somewhat Unimportant",
                                                                 "Somewhat Important",
                                                                 "Very Important")),
                                job  = factor(job, 
                                              levels = 1:4, 
                                              labels = c("Not Seeking Work",
                                                         "Looking for Work",
                                                         "Student",
                                                         "Employed"))
                                
)


## Create long dataset - categorical
ven.merge.cat <- ven.merge %>% 
  dplyr::select(lapop, age, marriage, education, race, salary, religion2, religiosity, job) %>%
  pivot_longer(cols = -lapop, names_to = "variable", values_to = "value") %>%
  drop_na()

bal.tab.cat.ven <- ven.merge.cat %>%
  group_by(variable) %>%
  do({
    
    unique_vals <- unique(.$value)
    pval_out <- map_df(unique_vals, function(x){
      x_lapop <- sum(.$value[.$lapop == 1] == x)
      x_hpz <- sum(.$value[.$lapop == 0] == x)
      n_lapop <- sum(.$lapop)
      n_hpz <- nrow(.) - n_lapop
      ptest_out <- prop.test(x = c(x_lapop, x_hpz), n = c(n_lapop, n_hpz))
      data.frame(value = x, 
                 mean_lapop = ptest_out$estimate[1], 
                 #se_lapop = sqrt(ptest_out$estimate[1] * (1 - ptest_out$estimate[1]) / n_lapop),
                 mean_hpz = ptest_out$estimate[2], 
                 #se_hpz = sqrt(ptest_out$estimate[2] * (1 - ptest_out$estimate[2]) / n_hpz),
                 diff_means = ptest_out$estimate[1] - ptest_out$estimate[2],
                 pval_difference = ptest_out$p.value) %>%
        mutate(value = as.factor(value)) 
    }) %>%
      arrange(value)
    
    pval_out  
  })  %>%
  rename(`Variable` = variable, 
         `Category` = value, 
         `LAPOP Mean` = mean_lapop, 
         #`LAPOP SE` = se_lapop, 
         `Sample Mean` = mean_hpz,
         #`Sample SE` = se_hpz,
         `Diff. in Means` = diff_means,
         `P-Value` = pval_difference) %>%
   mutate(Variable = case_when(Variable == "age"~"Age",
                              Variable == "education"~"Education",
                              Variable == "job"~"Employment",
                              Variable == "marriage"~"Marriage",
                              Variable == "race"~"Race",
                              Variable == "religion2"~"Religion",
                              Variable == "religiosity"~"Religiosity",
                              Variable == "salary"~"Salary"                              
                              )) %>%
  kable(align = "lccc",
        caption = "Comparing LAPOP and our Sample in Venezuela: Categorical Variables",
        label ="tab:lapop_ven_cat") %>%
  kable_styling(bootstrap_options = "striped") 


# show tables
bal.tab.num.ven %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 9)

bal.tab.cat.ven %>% 
  kable_styling(latex_options = c("HOLD_position"), font_size = 9)


@


\newpage
\section{Exploring what being left means for Venezuelans}
\label{SIsec:venleft}

For our Venezuelan respondents, this figure shows the relationship between being leftist and beliefs about privatization and support for Petro (if Venezuelans in Colombia could vote in the Colombian elections).
Support for Petro was not correlated with supporting the left, and identifying with the left is not correlated with the belief in government ownership of industries. 

<<ven_left_beliefs, eval = TRUE, echo = FALSE, tidy=TRUE, fig.width = 7, fig.height = 2.5, out.width= "1\\linewidth", fig.align='center', warning=FALSE, message=FALSE, strip.white=TRUE, fig.cap="For Venezuelans, association between being leftist and support for government ownership of industries (left plot), and supporting Petro (right plot).">>=

## Privatization
left_privatization_bi <- lm_robust(privatization ~ left_partisanship, data = df_ven)
# left_unemployed_ctrls <- lm_robust(update(pol_left ~ pol_guerilla, 
#                                     reformulate(c(".", labels(demo_ctrls_ven)))), data = df_col)

left_privatization_bi_plot <- left_privatization_bi %>%
        tidy() %>%
        filter(grepl("left_partisanship", term)) %>%
  ggplot(aes(x = term, y = estimate)) + 
  geom_point(size = 2, colour = "#D55E00") + 
  geom_label(aes(label = round(estimate, 2)), nudge_x = 0.2, colour = "#D55E00", size = 3) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), 
                width = 0, size = 1, colour = "#D55E00") + 
  scale_colour_manual(values = c("black")) + 
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  scale_y_continuous(breaks = c(-.5,0,.5), limits = c(-.6, .6)) + 
  scale_x_discrete(labels=str_wrap("Left ideology (1/0)", width = 10)) +
  labs(x = "", 
       y = "Estimate") +
  ggtitle(str_wrap("Government should own all industries (1-7)", width = 30)) +
  coord_flip() +
  yy_theme()

df_ven<-df_ven %>%
  mutate(petro_vote = case_when(pres_vote_col == 2~1,
                                pres_vote_col %in% c(0, 1, 3)~0,
                                TRUE~NA_real_)
         )

## Petro
left_petro_vote_bi <- lm_robust(petro_vote ~ left_partisanship, data = df_ven)
# left_unemployed_ctrls <- lm_robust(update(pol_left ~ pol_guerilla, 
#                                     reformulate(c(".", labels(demo_ctrls_ven)))), data = df_col)

left_petro_vote_bi_plot <- left_petro_vote_bi %>%
        tidy() %>%
        filter(grepl("left_partisanship", term)) %>%
  ggplot(aes(x = term, y = estimate)) + 
  geom_point(size = 2, colour = "#D55E00") + 
  geom_label(aes(label = round(estimate, 2)), nudge_x = 0.2, colour = "#D55E00", size = 3) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), 
                width = 0, size = 1, colour = "#D55E00") + 
  scale_colour_manual(values = c("black")) + 
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  scale_y_continuous(breaks = c(-.5,0,.5), limits = c(-.6, .6)) + 
  scale_x_discrete(labels=str_wrap("Left ideology (1/0)", width = 10)) +
  labs(x = "", 
       y = "Estimate") +
  ggtitle(str_wrap("Would vote for Petro (1/0)", width = 30)) +
  coord_flip() +
  yy_theme()

#plot
(left_privatization_bi_plot + left_petro_vote_bi_plot)
@


%\newpage
\section{Additional Attitudes of Colombians about Venezuelans}
\label{SIsec:colatts}

More than half of our Colombian respondents (\Sexpr{mean(df_col$labor_compete1==1, na.rm=T)*100}\%) said that they are competing with Venezuelans for work and
\Sexpr{mean(df_col$fisc_service==0, na.rm=T)*100}\% thought that it had become harder to obtain public services. Unlike poorer developing countries with large migrant flows, Colombia is a net fiscal contributor to hosting Venezuelan migrants and Colombians feel the strain:
\Sexpr{mean(df_col$fisc_aid==0, na.rm=T)*100}\% believed that the international community is not providing enough aid and 
\Sexpr{mean(df_col$fisc_tax==0, na.rm=T)*100}\% thought their own taxes will go up because of the migrants. 
Finally, despite cultural similarities, Colombians had concerns about their societal impact: only 
\Sexpr{mean(df_col$socio_integrate==1, na.rm=T)*100}\% said that Venezuelans have integrated successfully.


\newpage
\section{Correlation between Misperceptions and Welfare Concerns}
\label{SIsec:miswel}

For Colombian respondents, there is no relationship between believing the majority of Venezuelans are leftist and concerns over taxes increasing or concerns over access to government services, due to the presence of Venezuelans. 

<<regress_taxes, eval = TRUE, echo = FALSE, tidy=TRUE, fig.width = 7, fig.height = 2.5, out.width= "1\\linewidth", fig.align='left', warning=FALSE, message=FALSE, strip.white=TRUE, fig.cap="Correlations between political misperceptions (that the majority of Venezuelans are leftist) and concerns over taxes increasing or less access to government services.">>=

df_col <- df_col %>% # flip these two vars
 mutate(fisc_tax_0 = case_when(fisc_tax ==1~0, 
                  fisc_tax ==0~1, TRUE~NA_real_),
        fisc_service_0 = case_when(fisc_service ==1~0, 
                  fisc_service ==0~1, TRUE~NA_real_)
     ) 


# Ven support left - taxes go up
pol_left_fisc_tax_bi <- lm_robust(pol_left ~ fisc_tax_0, data = df_col)
pol_left_fisc_tax_ctrls <- lm_robust(update(pol_left ~ fisc_tax_0, 
                                    reformulate(c(".", labels(demo_ctrls_col)))), data = df_col)

pol_left_fisc_tax_ctrls_plot <- pol_left_fisc_tax_ctrls %>%
        tidy() %>%
        filter(grepl("fisc_tax_0", term)) %>%
  ggplot(aes(x = term, y = estimate)) + 
  geom_point(size = 2, colour = "#660066") + 
  geom_label(aes(label = round(estimate, 2)), nudge_x = 0.2, colour = "#660066", size = 3) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), 
                width = 0, size = 1, colour = "#660066") + 
  scale_colour_manual(values = c("black")) + 
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  scale_y_continuous(breaks = c(-.5,0,.5), limits = c(-.6, .6)) + 
  scale_x_discrete(labels=str_wrap("With Venezuelans, my taxes will go up (1/0)", width = 13)) +
  labs(x = "", 
       y = "Estimate") +
  ggtitle(str_wrap("Most Venezuelans in Colombia support the left (1/0)", width = 30)) +
  coord_flip() +
  yy_theme()

# Ven support left - less access to govt service
pol_left_fisc_service_bi <- lm_robust(pol_left ~ fisc_service_0, data = df_col)
pol_left_fisc_service_ctrls <- lm_robust(update(pol_left ~ fisc_service_0, 
                                    reformulate(c(".", labels(demo_ctrls_col)))), data = df_col)

pol_left_fisc_service_ctrls_plot <- pol_left_fisc_service_ctrls %>%
        tidy() %>%
        filter(grepl("fisc_service_0", term)) %>%
  ggplot(aes(x = term, y = estimate)) + 
  geom_point(size = 2, colour = "#660066") + 
  geom_label(aes(label = round(estimate, 2)), nudge_x = 0.2, colour = "#660066", size = 3) +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), 
                width = 0, size = 1, colour = "#660066") + 
  scale_colour_manual(values = c("black")) + 
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  scale_y_continuous(breaks = c(-.5,0,.5), limits = c(-.6, .6)) + 
  scale_x_discrete(labels=str_wrap("With Venezuelans, difficult to access govt services (1/0)", width = 13)) +
  labs(x = "", 
       y = "Estimate") +
  ggtitle(str_wrap("Most Venezuelans in Colombia support the left (1/0)", width = 30)) +
  coord_flip() +
  yy_theme()

#plot
(pol_left_fisc_tax_ctrls_plot + pol_left_fisc_service_ctrls_plot) 

@



\newpage
\section{Additional Conjoint Results}
\label{SIsec:conj}

\subsection{Main AMCE for Venezuelan respondents}

<<main_conj_ven, eval = TRUE, echo = FALSE, tidy=TRUE, fig.width = 8, fig.height = 7, out.width= "1\\linewidth", fig.align='left', warning=FALSE, message=FALSE, strip.white=TRUE, fig.cap="AMCE conjoint estimates with 95\\% CIs for Venezuelan respondents.">>=


#### Colombians

## Full conjoint dataframe
df_cjoint_col <- df_col %>% 
 dplyr::select(matches("con[1-5]"), 
        starts_with("rd_"),
        partisanship,
        skilled_labor,
        salary,
        benefits_bi,
        indir_contact_index_2,
        city,
        race) %>% 
 dplyr::mutate(id = 1:n(),
        partisanship = as.factor(partisanship),
        skilled_labor = as.factor(skilled_labor),
        salary = as.factor(salary),
        benefits_bi = as.factor(benefits_bi),
        indir_contact_index_2 = as.factor(indir_contact_index_2),
        city = as.factor(city),
        race = as.factor(race))

## Clean responses
df_cjoint_response_col <- df_cjoint_col %>%
 dplyr::select(id, matches("con[1-5]")) %>%
 gather(key = round, value = profile_choice, -id) %>%
 mutate(round = parse_number(round))

## Clean attributes
df_cjoint_covs_col <- df_cjoint_col %>%
 dplyr::select(id, starts_with("rd_")) %>%
 gather(key = feature_round_profile, value = attribute, -id) %>%
 mutate(round = parse_number(feature_round_profile),
     profile = case_when(grepl("_a_", feature_round_profile)~"a",
               grepl("_b_", feature_round_profile)~"b",
               TRUE~"ruh-roh"),
     feature = gsub("rd_[1-5]_[a|b]_", "", feature_round_profile)) %>%
 dplyr::select(-feature_round_profile) %>%
 spread(key = feature, value = attribute)

## Bind responses and attributes together
df_cjoint_clean_col <- inner_join(df_cjoint_response_col, df_cjoint_covs_col,
               by = c("id", "round")) %>%
 left_join(df_cjoint_col %>% dplyr::select(id, 
                      partisanship,
                      skilled_labor,
                      salary,
                      benefits_bi,
                      indir_contact_index_2,
                      city,
                      race)) %>%
 mutate(chosen_profile = case_when(profile_choice == 1 & profile == "a"~1,
                  profile_choice == 2 & profile == "b"~1,
                  profile_choice == 1 & profile == "b"~0,
                  profile_choice == 2 & profile == "a"~0,
                  TRUE~NA_real_),
     Partisanship = as.factor(Partisanship), 
     `Skill Level` = as.factor(`Skill Level`), 
     `Identity of Migrant` = as.factor(`Identity of Migrant`),
     `Probability of Employment` = as.factor(`Probability of Employment`), 
     Race = as.factor(Race), 
     `Reason for Leaving` = as.factor(`Reason for Leaving`), 
     Gender = as.factor(Gender)
     ) %>%
 drop_na(-c(partisanship,
        skilled_labor,
        salary,
        benefits_bi,
        indir_contact_index_2,
        city,
        race)) %>%
 mutate(pair_id = paste0(id, round)) %>%
 group_by(pair_id) %>%
 filter(n() == 2) %>% ungroup()

## Run conjoint estimator
set_baseline <- list(Partisanship = "Left",
           `Skill Level` = "Low",
           `Identity of Migrant` = "Colombians displaced from around Cali",
           `Probability of Employment` = "Cannot work",
           Race = "Mestizo",
           `Reason for Leaving` = "Fear of crime in their area")

df_cjoint_clean_col$Partisanship <- fct_relevel(df_cjoint_clean_col$Partisanship, "Left", "Center", "Right")
df_cjoint_clean_col$`Skill Level` <- fct_relevel(df_cjoint_clean_col$`Skill Level`, "Low", "Medium", "High")

amce_out.col <- cjoint::amce(chosen_profile ~ Partisanship + `Identity of Migrant` + `Skill Level` +
    `Probability of Employment` +`Reason for Leaving` + Race + Gender,
   data = df_cjoint_clean_col, 
   respondent.id = "id", 
   baselines = set_baseline)

amce_out.col$estimates <- lapply(amce_out.col$estimates, function(x){
 est <- x[1,]
 stderr <- x[2,]
 ci_low <- est - stderr * qnorm(.975)
 ci_hi <- est + stderr * qnorm(.975)
 xout <- rbind(est, stderr, ci_low, ci_hi)
 colnames(xout) <- colnames(x)
 return(xout)
})

## Plot results
cbPalette <- c("#99098d", "#000000", "#D55E00", "#009E73", "#E69F00", "#999999", "#0072B2")

amce_out.col.plot <- cjoint_plot(amce_out.col, 
   xlim=c(-.27,.18), breaks=c(-.2, -.1, 0, .1), labels=c("-.2", "-.1", "0",".1"), 
   xlab="Change in Pr(Preferred Migrant)",
   #main="Colombian respondents",
   lwd = 1.5,
   colors = cbPalette,
   group.order = c("Partisanship", "Identity of Migrant", "Skill Level",
           "Probability of Employment", "Reason for Leaving", "Race", "Gender")) +
   theme_bw() + 
   theme(legend.position = "none")


#### Venezuelans

## Full conjoint dataframe
df_cjoint_ven <- df_ven %>% 
 dplyr::select(matches("con[1-5]"), 
        starts_with("rd_"),
        partisanship,
        skilled_labor,
        salary,
        benefits_bi,
        city) %>% 
 dplyr::mutate(id = 1:n(),
        partisanship = as.factor(partisanship),
        skilled_labor = as.factor(skilled_labor),
        salary = as.factor(salary),
        benefits_bi = as.factor(benefits_bi),
        city = as.factor(city))

## Clean responses
df_cjoint_response_ven <- df_cjoint_ven %>%
 dplyr::select(id, matches("con[1-5]")) %>%
 gather(key = round, value = profile_choice, -id) %>%
 mutate(round = parse_number(round))

## Clean attributes
df_cjoint_covs_ven <- df_cjoint_ven %>%
 dplyr::select(id, starts_with("rd_")) %>%
 gather(key = feature_round_profile, value = attribute, -id) %>%
 mutate(round = parse_number(feature_round_profile),
     profile = case_when(grepl("_a_", feature_round_profile)~"a",
               grepl("_b_", feature_round_profile)~"b",
               TRUE~"ruh-roh"),
     feature = gsub("rd_[1-5]_[a|b]_", "", feature_round_profile)) %>%
 dplyr::select(-feature_round_profile) %>%
 spread(key = feature, value = attribute)

## Bind responses and attributes together
df_cjoint_clean_ven <- inner_join(df_cjoint_response_ven, df_cjoint_covs_ven,
               by = c("id", "round")) %>%
 left_join(df_cjoint_ven %>% dplyr::select(id,
                      partisanship,
                      skilled_labor,
                      salary,
                      benefits_bi,
                      city)) %>%
 mutate(chosen_profile = case_when(profile_choice == 1 & profile == "a"~1,
                  profile_choice == 2 & profile == "b"~1,
                  profile_choice == 1 & profile == "b"~0,
                  profile_choice == 2 & profile == "a"~0,
                  TRUE~NA_real_),
     Partisanship = as.factor(Partisanship), 
     `Skill Level` = as.factor(`Skill Level`), 
     `Identity of Migrant` = as.factor(`Identity of Migrant`),
     `Probability of Employment` = as.factor(`Probability of Employment`), 
     Race = as.factor(Race), 
     `Reason for Leaving` = as.factor(`Reason for Leaving`), 
     Gender = as.factor(Gender)
     ) %>%
 drop_na(-c(partisanship,
        skilled_labor,
        salary,
        benefits_bi,
        city)) %>%
 mutate(pair_id = paste0(id, round)) %>%
 group_by(pair_id) %>%
 filter(n() == 2) %>% ungroup()
 #drop_na(profile_choice)

## Run conjoint estimator
set_baseline <- list(Partisanship = "Left",
           `Skill Level` = "Low",
           `Identity of Migrant` = "Colombians displaced from around Cali",
           `Probability of Employment` = "Cannot work",
           Race = "Mestizo",
           `Reason for Leaving` = "Fear of crime in their area;")

df_cjoint_clean_ven$Partisanship <- fct_relevel(df_cjoint_clean_ven$Partisanship, "Left", "Center", "Right")
df_cjoint_clean_ven$`Skill Level` <- fct_relevel(df_cjoint_clean_ven$`Skill Level`, "Low", "Medium", "High")

amce_out.ven <- cjoint::amce(chosen_profile ~ Partisanship + `Identity of Migrant` + `Skill Level` +
    `Probability of Employment` +`Reason for Leaving` + Race + Gender,
   data = df_cjoint_clean_ven, 
   respondent.id = "id", 
   baselines = set_baseline)

amce_out.ven$estimates <- lapply(amce_out.ven$estimates, function(x){
 est <- x[1,]
 stderr <- x[2,]
 ci_low <- est - stderr * qnorm(.975)
 ci_hi <- est + stderr * qnorm(.975)
 xout <- rbind(est, stderr, ci_low, ci_hi)
 colnames(xout) <- colnames(x)
 return(xout)
})

## Plot results
cbPalette <- c("#99098d", "#000000", "#D55E00", "#009E73", "#E69F00", "#999999", "#0072B2")

amce_out.ven.plot <- cjoint_plot(amce_out.ven, 
   xlim=c(-.27,.18), breaks=c(-.2, -.1, 0, .1), labels=c("-.2", "-.1", "0",".1"), 
   xlab="Change in Pr(Preferred Migrant)",
   main="Venezuelan respondents",
   group.order = c("Partisanship", "Identity of Migrant", "Skill Level",
           "Probability of Employment", "Reason for Leaving", "Race", "Gender"),
   colors = cbPalette) +
   theme_bw() + 
   theme(legend.position = "none")

## plot
#amce_out.col.plot + amce_out.ven.plot + plot_layout(ncol=2, widths=c(1,1))

amce_out.ven.plot

@


\newpage
\subsection{Conditional AMCE on respondent characteristics}

<<cond_conj_ven, eval = TRUE, echo = FALSE, tidy=TRUE, fig.width = 9, fig.height = 5, out.width= "1\\linewidth", fig.align='center', warning=FALSE, message=FALSE, strip.white=TRUE, fig.cap="AMCE estimates with 95\\% CIs for Venezuelan respondents, of migrant profile political polarization conditional on respondent political polarization and city.">>=

#### Respondent political polarization

#### Colombians
amce_out.col <- cjoint::amce(chosen_profile ~ Partisanship:partisanship + `Identity of Migrant` + `Skill Level` +
       `Probability of Employment` +`Reason for Leaving` + Race + Gender,
     data = df_cjoint_clean_col, 
     respondent.id = "id", respondent.varying = "partisanship",
     baselines = set_baseline)

## Plot results
amce_out.col <- cjoint_plot(amce_out.col) 
amce_out.col <- amce_out.col$data %>% 
  filter(grepl("Partisanship", var)) %>%
  filter(printvar %in% c("   Center", "   Right")) %>%
  mutate(facet = as.character(facet),
         facet2 = case_when(facet == "Unconditional"~"Unconditional",
                           facet == "Conditional on\npartisanship = 1"~"Conditional on Left",
                           facet == "Conditional on\npartisanship = 2"~"Conditional on Center",
                           facet == "Conditional on\npartisanship = 3"~"Conditional on Right",
                           TRUE~"ruh-roh"),
         facet3 = fct_relevel(facet2, "Unconditional", 
                             "Conditional on Left", 
                             "Conditional on Center",
                             "Conditional on Right"),
         printvar = fct_relevel(printvar, 
                                "   Right",
                                "   Center"))

amce_out_pol.col.plot <- ggplot(amce_out.col, aes(printvar, pe)) + 
  geom_point(size = 2.5, colour = "#D55E00") + 
  geom_errorbar(aes(ymin = lower, ymax = upper), 
                width = 0, size = .5, colour = "#D55E00") + 
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  ylim(-.15, .2) +
  facet_wrap(~facet3, ncol = 1) + 
  #ggtitle("Colombian respondents") +
  ggtitle("By Ideology") +
  ylab("Change in Pr(Preferred Migrant)") +
  theme_bw() + 
  theme(legend.position = "none",
        axis.title.y = element_blank()) +
  coord_flip()

#### Venezuelans
amce_out.ven <- cjoint::amce(chosen_profile ~ Partisanship:partisanship + `Identity of Migrant` + `Skill Level` +
       `Probability of Employment` +`Reason for Leaving` + Race + Gender,
     data = df_cjoint_clean_ven[!is.na(df_cjoint_clean_ven$partisanship),], 
     respondent.id = "id", respondent.varying = "partisanship",
     baselines = set_baseline)

## Plot results
amce_out.ven <- cjoint_plot(amce_out.ven) 
amce_out.ven <- amce_out.ven$data %>% 
  filter(grepl("Partisanship", var)) %>%
  filter(printvar %in% c("   Center", "   Right")) %>%
  mutate(facet = as.character(facet),
         facet2 = case_when(facet == "Unconditional"~"Unconditional",
                           facet == "Conditional on\npartisanship = 1"~"Conditional on Left",
                           facet == "Conditional on\npartisanship = 2"~"Conditional on Center",
                           facet == "Conditional on\npartisanship = 3"~"Conditional on Right",
                           TRUE~"ruh-roh"),
         facet3 = fct_relevel(facet2, "Unconditional", 
                             "Conditional on Left", 
                             "Conditional on Center",
                             "Conditional on Right"),
         printvar = fct_relevel(printvar, 
                                "   Right",
                                "   Center"))

amce_out_pol.ven.plot <- ggplot(amce_out.ven, aes(printvar, pe)) + 
  geom_point(size = 2.5, colour = "#D55E00") + 
  geom_errorbar(aes(ymin = lower, ymax = upper),
                width = 0, size = .5, colour = "#D55E00") + 
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  ylim(-.15, .2) +
  facet_wrap(~facet3, ncol = 1) + 
  #ggtitle("Venezuelan respondents") +
  ggtitle("By Ideology") +
  ylab("Change in Pr(Preferred Migrant)") +
  theme_bw() + 
  theme(legend.position = "none",
        axis.title.y = element_blank(),
        ) +
  coord_flip()

#### Respondent city

#### Colombians
amce_out.col <- cjoint::amce(chosen_profile ~ Partisanship:city + `Identity of Migrant` + `Skill Level` +
       `Probability of Employment` +`Reason for Leaving` + Race + Gender,
     data = df_cjoint_clean_col, 
     respondent.id = "id", respondent.varying = "city",
     baselines = set_baseline)

## Plot results
amce_out.col <- cjoint_plot(amce_out.col) 
amce_out.col <- amce_out.col$data %>% 
  filter(grepl("Partisanship", var)) %>%
  filter(printvar %in% c("   Center", "   Right")) %>%
  mutate(facet = as.character(facet),
         facet2 = case_when(facet == "Unconditional"~"Unconditional",
                           facet == "Conditional on\ncity = Cali"~"Conditional on Cali",
                           facet == "Conditional on\ncity = Cúcuta"~"Conditional on Cúcuta",
                           TRUE~"ruh-roh"),
         facet3 = fct_relevel(facet2, "Unconditional", 
                             "Conditional on Cali", 
                             "Conditional on Cúcuta"),
         printvar = fct_relevel(printvar, 
                                "   Right",
                                "   Center"))

amce_out_city.col.plot <- ggplot(amce_out.col, aes(printvar, pe)) + 
  geom_point(size = 2.5, colour = "#D55E00") + 
  geom_errorbar(aes(ymin = lower, ymax = upper), 
                width = 0, size = .5, colour = "#D55E00") + 
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  ylim(-.1, .15) +
  facet_wrap(~facet3, ncol = 1) + 
  ggtitle("By City") +
  #ggtitle("Colombian respondents") +
  ylab("Change in Pr(Preferred Migrant)") +
  theme_bw() + 
  theme(legend.position = "none",
        axis.title.y = element_blank()) +
  coord_flip()

#### Venezuelans
amce_out.ven <- cjoint::amce(chosen_profile ~ Partisanship:city + `Identity of Migrant` + `Skill Level` +
       `Probability of Employment` +`Reason for Leaving` + Race + Gender,
     data = df_cjoint_clean_ven, 
     respondent.id = "id", respondent.varying = "city",
     baselines = set_baseline)

## Plot results
amce_out.ven <- cjoint_plot(amce_out.ven) 
amce_out.ven <- amce_out.ven$data %>% 
  filter(grepl("Partisanship", var)) %>%
  filter(printvar %in% c("   Center", "   Right")) %>%
  mutate(facet = as.character(facet),
         facet2 = case_when(facet == "Unconditional"~"Unconditional",
                           facet == "Conditional on\ncity = Cali"~"Conditional on Cali",
                           facet == "Conditional on\ncity = Cúcuta"~"Conditional on Cúcuta",
                           TRUE~"ruh-roh"),
         facet3 = fct_relevel(facet2, "Unconditional", 
                             "Conditional on Cali", 
                             "Conditional on Cúcuta"),
         printvar = fct_relevel(printvar, 
                                "   Right",
                                "   Center"))

amce_out_city.ven.plot <- ggplot(amce_out.ven, aes(printvar, pe)) + 
  geom_point(size = 2.5, colour = "#D55E00") + 
  geom_errorbar(aes(ymin = lower, ymax = upper),
                width = 0, size = .5, colour = "#D55E00") + 
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  ylim(-.1, .15) +
  facet_wrap(~facet3, ncol = 1) + 
  #ggtitle("Venezuelan respondents") +
  ggtitle("By City") +
  ylab("Change in Pr(Preferred Migrant)") +
  theme_bw() + 
  theme(legend.position = "none",
        axis.title.y = element_blank(),
        ) +
  coord_flip()


## plot
amce_out_pol.ven.plot + amce_out_city.ven.plot + plot_layout(ncol=2, widths=c(1,1))

@

<<cond_conj_skill, eval = TRUE, echo = FALSE, tidy=TRUE, fig.width = 8, fig.height = 5.5, out.width= "1\\linewidth", fig.align='center', warning=FALSE, message=FALSE, strip.white=TRUE, fig.cap="AMCE of migrant profile skill level conditional on respondent skill">>=

#### Respondent skilled_labor

#### Colombians
amce_out.col <- cjoint::amce(chosen_profile ~ Partisanship + `Identity of Migrant` + 
                       `Skill Level`:skilled_labor + `Probability of Employment` +
                       `Reason for Leaving` + Race + Gender,
     data = df_cjoint_clean_col[!is.na(df_cjoint_clean_col$skilled_labor),], 
     respondent.id = "id", respondent.varying = "skilled_labor",
     baselines = set_baseline)

## Plot results
amce_out.col <- cjoint_plot(amce_out.col) 
amce_out.col <- amce_out.col$data %>% 
  filter(grepl("SkillLevel", var)) %>%
  filter(printvar %in% c("   Medium", "   High")) %>%
  mutate(facet = as.character(facet),
         facet2 = case_when(facet == "Unconditional"~"Unconditional",
                           facet == "Conditional on\nskilled_labor = 0"~"Conditional on Unemployed",
                           facet == "Conditional on\nskilled_labor = 1"~"Conditional on Low skilled",
                           facet == "Conditional on\nskilled_labor = 2"~"Conditional on High skilled",
                           TRUE~"ruh-roh"),
         facet3 = fct_relevel(facet2, "Unconditional", 
                             "Conditional on Unemployed", 
                             "Conditional on Low skilled",
                             "Conditional on High skilled"),
         printvar = fct_relevel(printvar, 
                                "   High",
                                "   Medium"))

amce_out_skill.col.plot <- ggplot(amce_out.col, aes(printvar, pe)) + 
  geom_point(size = 2.5, colour = "#0072B2") + 
  geom_errorbar(aes(ymin = lower, ymax = upper),
                width = 0, size = .5, colour = "#0072B2") + 
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  ylim(-.15, .2) +
  facet_wrap(~facet3, ncol = 1) + 
  ggtitle("Colombian respondents") +
  ylab("Change in Pr(Preferred Migrant)") +
  theme_bw() + 
  theme(legend.position = "none",
        axis.title.y = element_blank()) +
  coord_flip()

#### Venezuelans
amce_out.ven <- cjoint::amce(chosen_profile ~ Partisanship + `Identity of Migrant` + 
                       `Skill Level`:skilled_labor + `Probability of Employment` +
                       `Reason for Leaving` + Race + Gender,
     data = df_cjoint_clean_ven[!is.na(df_cjoint_clean_ven$skilled_labor),], 
     respondent.id = "id", respondent.varying = "skilled_labor",
     baselines = set_baseline)

## Plot results
amce_out.ven <- cjoint_plot(amce_out.ven) 
amce_out.ven <- amce_out.ven$data %>% 
  filter(grepl("SkillLevel", var)) %>%
  filter(printvar %in% c("   Medium", "   High")) %>%
  mutate(facet = as.character(facet),
         facet2 = case_when(facet == "Unconditional"~"Unconditional",
                           facet == "Conditional on\nskilled_labor = 0"~"Conditional on Unemployed",
                           facet == "Conditional on\nskilled_labor = 1"~"Conditional on Low skilled",
                           facet == "Conditional on\nskilled_labor = 2"~"Conditional on High skilled",
                           TRUE~"ruh-roh"),
         facet3 = fct_relevel(facet2, "Unconditional", 
                             "Conditional on Unemployed", 
                             "Conditional on Low skilled",
                             "Conditional on High skilled"),
         printvar = fct_relevel(printvar, 
                                "   High",
                                "   Medium"))

amce_out_skill.ven.plot <- ggplot(amce_out.ven, aes(printvar, pe)) + 
  geom_point(size = 2.5, colour = "#0072B2") + 
  geom_errorbar(aes(ymin = lower, ymax = upper),
                width = 0, size = .5, colour = "#0072B2") + 
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  ylim(-.15, .2) +
  facet_wrap(~facet3, ncol = 1) + 
  ggtitle("Venezuelan respondents") +
  ylab("Change in Pr(Preferred Migrant)") +
  theme_bw() + 
  theme(legend.position = "none",
        axis.title.y = element_blank(),
        axis.text.y = element_blank()
        ) +
  coord_flip()


## plot
amce_out_skill.col.plot + amce_out_skill.ven.plot + plot_layout(ncol=2, widths=c(1,1))

@

<<cond_conj_salary, eval = TRUE, echo = FALSE, tidy=TRUE, fig.width = 8, fig.height = 5.5, out.width= "1\\linewidth", fig.align='center', warning=FALSE, message=FALSE, strip.white=TRUE, fig.cap="AMCE of migrant profile ability to find work conditional on respondent salary.">>=

#### Respondent salary

#### Colombians
amce_out.col <- cjoint::amce(chosen_profile ~ Partisanship + `Identity of Migrant` + 
                       `Skill Level` + `Probability of Employment`:salary +
                       `Reason for Leaving` + Race + Gender,
     data = df_cjoint_clean_col[!is.na(df_cjoint_clean_col$salary),], 
     respondent.id = "id", respondent.varying = "salary",
     baselines = set_baseline)

## Plot results
amce_out.col <- cjoint_plot(amce_out.col) 
amce_out.col <- amce_out.col$data %>% 
  filter(grepl("ProbabilityofEmployment", var)) %>%
  filter(printvar %in% c("   Unlikely to find work", "   Very likely to find work")) %>%
  mutate(facet = as.character(facet),
         facet2 = case_when(facet == "Unconditional"~"Unconditional",
                           facet == "Conditional on\nsalary = 1"~"Conditional on Not Enough, Major Difficulties",
                           facet == "Conditional on\nsalary = 2"~"Conditional on Not Enough, Some Difficulties",
                           facet == "Conditional on\nsalary = 3"~"Conditional on Enough, Some Difficulties",
                           facet == "Conditional on\nsalary = 4"~"Conditional on Enough with Savings",
                           TRUE~"ruh-roh"),
         facet3 = fct_relevel(facet2, "Unconditional", 
                             "Conditional on Not Enough, Major Difficulties", 
                             "Conditional on Not Enough, Some Difficulties",
                             "Conditional on Enough, Some Difficulties",
                             "Conditional on Enough with Savings"),
         printvar = fct_relevel(printvar, 
                                "   Very likely to find work",
                                "   Unlikely to find work"))

amce_out_salary.col.plot <- ggplot(amce_out.col, aes(printvar, pe)) + 
  geom_point(size = 2.5, colour = "#009E73") + 
  geom_errorbar(aes(ymin = lower, ymax = upper), 
                width = 0, size = .5, colour = "#009E73") + 
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  ylim(-.15, .2) +
  facet_wrap(~facet3, ncol = 1) + 
  ggtitle("Colombian respondents") +
  ylab("Change in Pr(Preferred Migrant)") +
  theme_bw() + 
  theme(legend.position = "none",
        axis.title.y = element_blank()) +
  coord_flip()

#### Venezuelans
amce_out.ven <- cjoint::amce(chosen_profile ~ Partisanship + `Identity of Migrant` + 
                       `Skill Level` + `Probability of Employment`:salary +
                       `Reason for Leaving` + Race + Gender,
     data = df_cjoint_clean_ven[!is.na(df_cjoint_clean_ven$salary),], 
     respondent.id = "id", respondent.varying = "salary",
     baselines = set_baseline)

## Plot results
amce_out.ven <- cjoint_plot(amce_out.ven) 
amce_out.ven <- amce_out.ven$data %>% 
  filter(grepl("ProbabilityofEmployment", var)) %>%
  filter(printvar %in% c("   Unlikely to find work", "   Very likely to find work")) %>%
  mutate(facet = as.character(facet),
         facet2 = case_when(facet == "Unconditional"~"Unconditional",
                           facet == "Conditional on\nsalary = 1"~"Conditional on Not Enough, Major Difficulties",
                           facet == "Conditional on\nsalary = 2"~"Conditional on Not Enough, Some Difficulties",
                           facet == "Conditional on\nsalary = 3"~"Conditional on Enough, Some Difficulties",
                           facet == "Conditional on\nsalary = 4"~"Conditional on Enough with Savings",
                           TRUE~"ruh-roh"),
         facet3 = fct_relevel(facet2, "Unconditional", 
                             "Conditional on Not Enough, Major Difficulties", 
                             "Conditional on Not Enough, Some Difficulties",
                             "Conditional on Enough, Some Difficulties",
                             "Conditional on Enough with Savings"),
         printvar = fct_relevel(printvar, 
                                "   Very likely to find work",
                                "   Unlikely to find work"))

amce_out_salary.ven.plot <- ggplot(amce_out.ven, aes(printvar, pe)) + 
  geom_point(size = 2.5, colour = "#009E73") + 
  geom_errorbar(aes(ymin = lower, ymax = upper), 
                width = 0, size = .5, colour = "#009E73") + 
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  ylim(-.15, .2) +
  facet_wrap(~facet3, ncol = 1) + 
  ggtitle("Venezuelan respondents") +
  ylab("Change in Pr(Preferred Migrant)") +
  theme_bw() + 
  theme(legend.position = "none",
        axis.title.y = element_blank(),
        axis.text.y = element_blank()
        ) +
  coord_flip()


## plot
amce_out_salary.col.plot + amce_out_salary.ven.plot + plot_layout(ncol=2, widths=c(1,1))


@

<<cond_conj_benefits, eval = TRUE, echo = FALSE, tidy=TRUE, fig.width = 8, fig.height = 4, out.width= "1\\linewidth", fig.align='center', warning=FALSE, message=FALSE, strip.white=TRUE, fig.cap="AMCE of migrant profile ability to find work conditional on respondent public benefits.">>=

#### Respondent benefits

#### Colombians
amce_out.col <- cjoint::amce(chosen_profile ~ Partisanship + `Identity of Migrant` + 
                       `Skill Level` + `Probability of Employment`:benefits_bi +
                       `Reason for Leaving` + Race + Gender,
     data = df_cjoint_clean_col[!is.na(df_cjoint_clean_col$benefits_bi),], 
     respondent.id = "id", respondent.varying = "benefits_bi",
     baselines = set_baseline)

## Plot results
amce_out.col <- cjoint_plot(amce_out.col) 
amce_out.col <- amce_out.col$data %>% 
  filter(grepl("ProbabilityofEmployment", var)) %>%
  filter(printvar %in% c("   Unlikely to find work", "   Very likely to find work")) %>%
  mutate(facet = as.character(facet),
         facet2 = case_when(facet == "Unconditional"~"Unconditional",
                           facet == "Conditional on\nbenefits_bi = 0"~"Conditional on Not Receiving Public Benefits",
                           facet == "Conditional on\nbenefits_bi = 1"~"Conditional on Receiving Public Benefits",
                           TRUE~"ruh-roh"),
         facet3 = fct_relevel(facet2, "Unconditional", 
                             "Conditional on Not Receiving Public Benefits", 
                             "Conditional on Receiving Public Benefits"),
         printvar = fct_relevel(printvar, 
                                "   Very likely to find work",
                                "   Unlikely to find work"))

amce_out_benefits.col.plot <- ggplot(amce_out.col, aes(printvar, pe)) + 
  geom_point(size = 2.5, colour = "#009E73") + 
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0, 
                size = .5, colour = "#009E73") + 
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  ylim(-.15, .2) +
  facet_wrap(~facet3, ncol = 1) + 
  ggtitle("Colombian respondents") +
  ylab("Change in Pr(Preferred Migrant)") +
  theme_bw() + 
  theme(legend.position = "none",
        axis.title.y = element_blank()) +
  coord_flip()

#### Venezuelans
amce_out.ven <- cjoint::amce(chosen_profile ~ Partisanship + `Identity of Migrant` + 
                       `Skill Level` + `Probability of Employment`:benefits_bi +
                       `Reason for Leaving` + Race + Gender,
     data = df_cjoint_clean_ven[!is.na(df_cjoint_clean_ven$benefits_bi),], 
     respondent.id = "id", respondent.varying = "benefits_bi",
     baselines = set_baseline)

## Plot results
amce_out.ven <- cjoint_plot(amce_out.ven) 
amce_out.ven <- amce_out.ven$data %>% 
  filter(grepl("ProbabilityofEmployment", var)) %>%
  filter(printvar %in% c("   Unlikely to find work", "   Very likely to find work")) %>%
  mutate(facet = as.character(facet),
         facet2 = case_when(facet == "Unconditional"~"Unconditional",
                           facet == "Conditional on\nbenefits_bi = 0"~"Conditional on Not Receiving Public Benefits",
                           facet == "Conditional on\nbenefits_bi = 1"~"Conditional on Receiving Public Benefits",
                           TRUE~"ruh-roh"),
         facet3 = fct_relevel(facet2, "Unconditional", 
                             "Conditional on Not Receiving Public Benefits", 
                             "Conditional on Receiving Public Benefits"),
         printvar = fct_relevel(printvar, 
                                "   Very likely to find work",
                                "   Unlikely to find work"))

amce_out_benefits.ven.plot <- ggplot(amce_out.ven, aes(printvar, pe)) + 
  geom_point(size = 2.5, colour = "#009E73") + 
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0, 
                size = .5, colour = "#009E73") + 
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  ylim(-.15, .2) +
  facet_wrap(~facet3, ncol = 1) + 
  ggtitle("Venezuelan respondents") +
  ylab("Change in Pr(Preferred Migrant)") +
  theme_bw() + 
  theme(legend.position = "none",
        axis.title.y = element_blank(),
        axis.text.y = element_blank()
        ) +
  coord_flip()


## plot
amce_out_benefits.col.plot + amce_out_benefits.ven.plot + plot_layout(ncol=2, widths=c(1,1))


@

<<cond_conj_indir_contact, eval = TRUE, echo = FALSE, tidy=TRUE, fig.width = 8, fig.height = 5, out.width= "1\\linewidth", fig.align='center', warning=FALSE, message=FALSE, strip.white=TRUE, fig.cap="AMCE of migrant profile ability to find work conditional on respondent frequency of indirect contact with Venezuelans (Colombians respondents only).">>=

#### Respondent Indirect Contact

#### Colombians
amce_out.col <- cjoint::amce(chosen_profile ~ Partisanship + `Identity of Migrant` + 
                       `Skill Level` + `Probability of Employment`:indir_contact_index_2 +
                       `Reason for Leaving` + Race + Gender,
     data = df_cjoint_clean_col[!is.na(df_cjoint_clean_col$indir_contact_index_2),], 
     respondent.id = "id", respondent.varying = "indir_contact_index_2",
     baselines = set_baseline)

## Plot results
amce_out.col <- cjoint_plot(amce_out.col) 
amce_out.col <- amce_out.col$data %>% 
  filter(grepl("ProbabilityofEmployment", var)) %>%
  filter(printvar %in% c("   Unlikely to find work", "   Very likely to find work")) %>%
  mutate(facet = as.character(facet),
         facet2 = case_when(facet == "Unconditional"~"Unconditional",
                           facet == "Conditional on\nindir_contact_index_2 = 0"~"Conditional on No Contact",
                           facet == "Conditional on\nindir_contact_index_2 = 1"~"Conditional on Some Contact",
                           facet == "Conditional on\nindir_contact_index_2 = 2"~"Conditional on Frequent Contact",
                           facet == "Conditional on\nindir_contact_index_2 = 3"~"Conditional on Contact All the Time",
                           TRUE~"ruh-roh"),
         facet3 = fct_relevel(facet2, "Unconditional", 
                             "Conditional on No Contact", 
                             "Conditional on Some Contact",
                             "Conditional on Frequent Contact",
                             "Conditional on Contact All the Time"),
         printvar = fct_relevel(printvar, 
                               "   Very likely to find work",
                               "   Unlikely to find work"))

amce_out_benefits.col.plot <- ggplot(amce_out.col, aes(printvar, pe)) + 
  geom_point(size = 2.5, colour = "#009E73") + 
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0, 
                size = .5, colour = "#009E73") + 
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  #ylim(-.15, .5) +
  facet_wrap(~facet3, ncol = 1) + 
  ggtitle("Colombian respondents") +
  ylab("Change in Pr(Preferred Migrant)") +
  theme_bw() + 
  theme(legend.position = "none",
        axis.title.y = element_blank()) +
  coord_flip()


## plot
amce_out_benefits.col.plot + plot_spacer() + plot_layout(ncol=2, widths=c(1,1))

@

<<cond_conj_origin_city, eval = TRUE, echo = FALSE, tidy=TRUE, fig.width = 8, fig.height = 4, out.width= "1\\linewidth", fig.align='center', warning=FALSE, message=FALSE, strip.white=TRUE, fig.cap="AMCE of migrant profile origin conditions conditional on respondent region (Colombians respondents only).">>=

#### Respondent benefits

#### Colombians
amce_out.col <- cjoint::amce(chosen_profile ~ Partisanship + `Identity of Migrant`:city + 
                       `Skill Level` + `Probability of Employment` +
                       `Reason for Leaving` + Race + Gender,
     data = df_cjoint_clean_col[!is.na(df_cjoint_clean_col$city),], 
     respondent.id = "id", respondent.varying = "city",
     baselines = set_baseline)

## Plot results
amce_out.col <- cjoint_plot(amce_out.col) 
amce_out.col <- amce_out.col$data %>% 
  filter(grepl("IdentityofMigrant", var)) %>%
  filter(printvar %in% c("   Colombians displaced from around Cúcuta", "   Venezuelans")) %>%
  mutate(facet = as.character(facet),
         facet2 = case_when(facet == "Unconditional"~"Unconditional",
                           facet == "Conditional on\ncity = Cali"~"Conditional on Cali",
                           facet == "Conditional on\ncity = Cúcuta"~"Conditional on Cúcuta",
                           TRUE~"ruh-roh"),
         facet3 = fct_relevel(facet2, "Unconditional", 
                             "Conditional on Cali", 
                             "Conditional on Cúcuta"),
         printvar = fct_relevel(printvar, 
                                "   Venezuelans",
                                "   Colombians displaced from around Cúcuta"))

amce_out_city.col.plot <- ggplot(amce_out.col, aes(printvar, pe)) + 
  geom_point(size = 2.5, colour = "#000000") + 
  geom_errorbar(aes(ymin = lower, ymax = upper),width = 0, 
                size = .5, colour = "#000000") + 
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  #ylim(-.15, .5) +
  facet_wrap(~facet3, ncol = 1) + 
  ggtitle("Colombian respondents") +
  ylab("Change in Pr(Preferred Migrant)") +
  theme_bw() + 
  theme(legend.position = "none",
        axis.title.y = element_blank()) +
  coord_flip()


## plot
amce_out_city.col.plot + plot_spacer() + plot_layout(ncol=2, widths=c(1,1))

@

%%% Race, conditioning on race, for Colombians only
<<cond_conj_race, eval = TRUE, echo = FALSE, tidy=TRUE, fig.width = 9, fig.height = 6, out.width= "1\\linewidth", fig.align='center', warning=FALSE, message=FALSE, strip.white=TRUE, fig.cap="AMCE estimates with 95\\% CIs of migrant profiles, conditional on respondent's race (Colombians respondents only).">>=

#### Colombians
amce_out.col <- cjoint::amce(chosen_profile ~ Partisanship:race + `Identity of Migrant`:race + 
                       `Skill Level`:race + `Probability of Employment`:race +
                       `Reason for Leaving`:race + Race:race + Gender:race,
     data = df_cjoint_clean_col[!is.na(df_cjoint_clean_col$race),], 
     respondent.id = "id", respondent.varying = "race",
     baselines = set_baseline)

## Plot results
amce_out.col <- cjoint_plot(amce_out.col, group.order = c("Partisanship", "Identity of Migrant", "Skill Level",
           "Probability of Employment", "Reason for Leaving", "Race", "Gender")) 

x_levels <- layer_scales(amce_out.col, 1, 1)$x$labels

amce_out.col <- amce_out.col$data %>% 
  #filter(grepl("Race", var)) %>%
  #filter(printvar %in% c("   (Baseline = Mestizo)", "   African descent")) %>%
  filter(printvar %in% x_levels) %>%
  mutate(facet = as.character(facet),
         facet2 = case_when(facet == "Unconditional"~"Unconditional",
                           facet == "Conditional on\nrace = 1"~"Conditional on White",
                           facet == "Conditional on\nrace = 2"~"Conditional on Mestizo",
                           facet == "Conditional on\nrace = 3"~"Conditional on Indigenous",
                           facet == "Conditional on\nrace = 4"~"Conditional on Afro/Black",
                           facet == "Conditional on\nrace = 5"~"Conditional on Mulato",
                           facet == "Conditional on\nrace = 6"~"Conditional on Other",
                           TRUE~"ruh-roh"),
         facet3 = fct_relevel(facet2, "Unconditional", 
                             "Conditional on White", 
                             "Conditional on Mestizo",
                             "Conditional on Afro/Black"),
         printvar = fct_relevel(printvar, x_levels)) %>%
  filter(!(facet3 %in% c("Conditional on Indigenous", "Conditional on Mulato", "Conditional on Other")))
  

amce_out_race.col.plot <- ggplot(amce_out.col, aes(printvar, pe, group = group, color = group)) + 
  geom_point() + 
  geom_errorbar(aes(ymin = lower, ymax = upper, color = group),
                width = 0, size = 1) + 
  scale_color_manual(values = cbPalette) +
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  ylim(-.35, .35) +
  facet_wrap(~facet3, ncol = 5) + 
  ggtitle("Colombian respondents") +
  ylab("Change in Pr(Migrant Preferred to Stay\nin the City and receive Benefits)") +
  theme_bw() + 
  theme(legend.position = "none",
        axis.title.y = element_blank()) +
  coord_flip()

amce_out_race.col.plot

@


%%% Political ideo, conditioning on Venezuelan friends and begging, for Colombians only
<<cond_conj_ven_friends, eval = FALSE, echo = FALSE, tidy=TRUE, fig.width = 7, fig.height = 5, out.width= ".8\\linewidth", fig.align='center', warning=FALSE, message=FALSE, strip.white=TRUE, fig.cap="AMCE estimates with 95\\% CIs of migrant profile political polarization, conditional on respondent's number of Venezuelan friends (Colombians respondents only).">>=

#### Respondent Venezuelan friends
amce_out.col <- cjoint::amce(chosen_profile ~ Partisanship:ven_friends + `Identity of Migrant` + `Skill Level` +
       `Probability of Employment` +`Reason for Leaving` + Race + Gender,
     data = df_cjoint_clean_col[!is.na(df_cjoint_clean_col$ven_friends),], 
     respondent.id = "id", respondent.varying = "ven_friends",
     baselines = set_baseline)

## Plot results
amce_out.col <- cjoint_plot(amce_out.col) 
amce_out.col <- amce_out.col$data %>% 
  filter(grepl("Partisanship", var)) %>%
  filter(printvar %in% c("   Center", "   Right")) %>%
  mutate(facet = as.character(facet),
         facet2 = case_when(facet == "Unconditional"~"Unconditional",
                           facet == "Conditional on\nven_friends = 0"~"Conditional on 0 friends",
                           facet == "Conditional on\nven_friends = 1"~"Conditional on 1-2 friends",
                           facet == "Conditional on\nven_friends = 2"~"Conditional on 3-5 friends",
                           facet == "Conditional on\nven_friends = 3"~"Conditional on 5+ friends",
                           TRUE~"ruh-roh"),
         facet3 = fct_relevel(facet2, "Unconditional", 
                             "Conditional on 0 friends",
                             "Conditional on 1-2 friends", 
                             "Conditional on 3-5 friends",
                             "Conditional on 5+ friends"),
         printvar = fct_relevel(printvar, 
                                "   Right",
                                "   Center"))

amce_out_ven_friends.col.plot <- ggplot(amce_out.col, aes(printvar, pe)) + 
  geom_point(size = 2.5, colour = "#D55E00") + 
  geom_errorbar(aes(ymin = lower, ymax = upper), 
                width = 0, size = .5, colour = "#D55E00") + 
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  ylim(-.1, .2) +
  facet_wrap(~facet3, ncol = 1) + 
  ggtitle("Colombian respondents") +
  ylab("Change in Pr(Preferred Migrant)") +
  theme_bw() + 
  theme(legend.position = "none",
        axis.title.y = element_blank()) +
  coord_flip()

#### Respondent Venezuelan begging
amce_out.col <- cjoint::amce(chosen_profile ~ Partisanship:ven_beg_bi + `Identity of Migrant` + `Skill Level` +
       `Probability of Employment` +`Reason for Leaving` + Race + Gender,
     data = df_cjoint_clean_col[!is.na(df_cjoint_clean_col$ven_beg_bi),], 
     respondent.id = "id", respondent.varying = "ven_beg_bi",
     baselines = set_baseline)

## Plot results
amce_out.col <- cjoint_plot(amce_out.col) 
amce_out.col <- amce_out.col$data %>% 
  filter(grepl("Partisanship", var)) %>%
  filter(printvar %in% c("   Center", "   Right")) %>%
  mutate(facet = as.character(facet),
         facet2 = case_when(facet == "Unconditional"~"Unconditional",
                           facet == "Conditional on\nven_beg_bi = 0"~"Conditional on not seeing\n Venezuelans begging every day",
                           facet == "Conditional on\nven_beg_bi = 1"~"Conditional on seeing\n Venezuelans begging every day",
                           TRUE~"ruh-roh"),
         facet3 = fct_relevel(facet2, "Unconditional", 
                             "Conditional on not seeing\n Venezuelans begging every day",
                             "Conditional on seeing\n Venezuelans begging every day"),
         printvar = fct_relevel(printvar, 
                                "   Right",
                                "   Center"))

amce_out_ven_beg.col.plot <- ggplot(amce_out.col, aes(printvar, pe)) + 
  geom_point(size = 2.5, colour = "#D55E00") + 
  geom_errorbar(aes(ymin = lower, ymax = upper), 
                width = 0, size = .5, colour = "#D55E00") + 
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  ylim(-.1, .3) +
  facet_wrap(~facet3, ncol = 1) + 
  ggtitle("Colombian respondents") +
  ylab("Change in Pr(Preferred Migrant)") +
  theme_bw() + 
  theme(legend.position = "none",
        axis.title.y = element_blank()) +
  coord_flip()

## plot
amce_out_ven_friends.col.plot + amce_out_ven_beg.col.plot + plot_layout(ncol=2, widths=c(1,1))

@

\newpage
\subsection{Interaction AMCE with profile attributes}
\label{SIsec:acie}

<<int_conj_part_reason, eval = FALSE, echo = FALSE, tidy=TRUE, fig.width = 8, fig.height = 3.5, out.width= "1\\linewidth", fig.align='center', warning=FALSE, message=FALSE, strip.white=TRUE, fig.cap="Interaction AMCE between partisanship and reason for leaving.">>=

cbPalette <- c("#99098d", "#000000", "#D55E00", "#009E73", "#E69F00", "#999999", "#0072B2")

#### Colombians

## Run interactive model
cjoint_inter_col <- CausalANOVA(chosen_profile ~ Partisanship + `Identity of Migrant` + 
                       `Skill Level` + `Probability of Employment` +
                       `Reason for Leaving` + Race + Gender,
       int2.formula = ~Partisanship:`Identity of Migrant` + Partisanship:`Reason for Leaving`,
       data = df_cjoint_clean_col, pair.id = df_cjoint_clean_col$pair_id, diff = TRUE,
       cluster = as.character(df_cjoint_clean_col$id), nway = 2, verbose = FALSE)

## Effect of partisanship conditional on reason for leaving
ceff_pol_reason_col <- ConditionalEffect(cjoint_inter_col, 
                                     treat.fac = "Partisanship", 
                                     cond.fac = "`Reason for Leaving`", 
                                     base.ind = 1, verbose = FALSE)

ceff_pol_reason_col <- bind_rows(
  ceff_pol_reason_col$ConditionalEffects %>% .[[1]] %>%  
    as.data.frame() %>%
    mutate(`Reason for Leaving` = "Fear of arrest by the government", Partisanship = rownames(.)),
  ceff_pol_reason_col$ConditionalEffects %>% .[[2]] %>%  
    as.data.frame() %>%
    mutate(`Reason for Leaving` = "Fear of crime in their area", Partisanship = rownames(.)),
  ceff_pol_reason_col$ConditionalEffects %>% .[[3]] %>%  
    as.data.frame() %>%
    mutate(`Reason for Leaving` = "Poverty", Partisanship = rownames(.)), 
  ceff_pol_reason_col$ConditionalEffects %>% .[[4]] %>%  
    as.data.frame() %>%
    mutate(`Reason for Leaving` = "Violence by guerillas, irregular forces, paramilitaries", 
           Partisanship = rownames(.)), 
) %>% filter(Partisanship != "Left")

ceff_pol_reason_col$Partisanship <- fct_relevel(ceff_pol_reason_col$Partisanship, "Right", "Center")

amcie_out_partisanship_reason_col.plot <- ggplot(ceff_pol_reason_col, 
                                            aes(Partisanship, 
                                                ConditionalEffect, 
                                                group = `Reason for Leaving`, 
                                                colour = `Reason for Leaving`)) + 
  geom_point(size = 2.5, position = position_dodge(.4)) + 
  geom_errorbar(aes(ymin = `2.5% CI`, ymax = `97.5% CI`), width = 0, 
                size = .5, position = position_dodge(.4)) + 
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  ylim(-.05, .2) +
  ggtitle("Colombian respondents") +
  ylab("Change in Pr(Preferred Migrant)") +
  theme_bw() +
  scale_colour_manual(values=cbPalette) +
  theme(legend.position="none",
        axis.title.y = element_blank()
        ) +
  coord_flip()


#### Venezuelans
## Run interactive model
cjoint_inter_ven <- CausalANOVA(chosen_profile ~ Partisanship + `Identity of Migrant` + 
                       `Skill Level` + `Probability of Employment` +
                       `Reason for Leaving` + Race + Gender,
       int2.formula = ~Partisanship:`Identity of Migrant` + Partisanship:`Reason for Leaving`,
       data = df_cjoint_clean_ven, pair.id = df_cjoint_clean_ven$pair_id, diff = TRUE,
       cluster = as.character(df_cjoint_clean_ven$id), nway = 2, verbose = FALSE)

## Effect of partisanship conditional on reason for leaving
ceff_pol_reason_ven <- ConditionalEffect(cjoint_inter_ven, 
                                     treat.fac = "Partisanship", 
                                     cond.fac = "`Reason for Leaving`", 
                                     base.ind = 1, verbose = FALSE)

ceff_pol_reason_ven <- bind_rows(
  ceff_pol_reason_ven$ConditionalEffects %>% .[[1]] %>%  
    as.data.frame() %>%
    mutate(`Reason for Leaving` = "Fear of arrest by the government", Partisanship = rownames(.)),
  ceff_pol_reason_ven$ConditionalEffects %>% .[[2]] %>%  
    as.data.frame() %>%
    mutate(`Reason for Leaving` = "Fear of crime in their area", Partisanship = rownames(.)),
  ceff_pol_reason_ven$ConditionalEffects %>% .[[3]] %>%  
    as.data.frame() %>%
    mutate(`Reason for Leaving` = "Poverty", Partisanship = rownames(.)), 
  ceff_pol_reason_ven$ConditionalEffects %>% .[[4]] %>%  
    as.data.frame() %>%
    mutate(`Reason for Leaving` = "Violence by guerillas, irregular forces, paramilitaries", 
           Partisanship = rownames(.)), 
) %>% filter(Partisanship != "Left")

ceff_pol_reason_ven$`Reason for Leaving`<-sub("by", "by\n", ceff_pol_reason_ven$`Reason for Leaving`)  
ceff_pol_reason_ven$`Reason for Leaving`<-sub("crime", "crime\n", ceff_pol_reason_ven$`Reason for Leaving`)  
ceff_pol_reason_ven$`Reason for Leaving`<-sub("guerrillas, ", "guerrillas,\n", ceff_pol_reason_ven$`Reason for Leaving`)  
ceff_pol_reason_ven$`Reason for Leaving`<-sub("irregular", "irregular\n",
                                             ceff_pol_reason_ven$`Reason for Leaving`)  
ceff_pol_reason_ven$`Reason for Leaving`<-sub("forces,", "forces,\n",
                                             ceff_pol_reason_ven$`Reason for Leaving`)  

ceff_pol_reason_ven$Partisanship <- fct_relevel(ceff_pol_reason_ven$Partisanship, "Right", "Center")

amcie_out_partisanship_reason_ven.plot <- ggplot(ceff_pol_reason_ven, 
                                            aes(Partisanship, 
                                                ConditionalEffect, 
                                                group = `Reason for Leaving`, 
                                                colour = `Reason for Leaving`)) + 
  geom_point(size = 2.5, position = position_dodge(.4)) + 
  geom_errorbar(aes(ymin = `2.5% CI`, ymax = `97.5% CI`), width = 0, 
                size = .5, position = position_dodge(.4)) + 
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  ylim(-.05, .2) +
  ggtitle("Venezuelan respondents") +
  ylab("Change in Pr(Preferred Migrant)") +
  theme_bw() +
  scale_colour_manual(values=cbPalette) +
  theme_bw() +
  theme(legend.justification=c(1,0), 
        legend.position=c(.99,0.01),
        legend.title = element_blank(),
        legend.key.size = unit(0.1, "cm"),
        axis.title.y = element_blank(),
        axis.text.y = element_blank()) +
  coord_flip()


## plot
amcie_out_partisanship_reason_col.plot + 
  amcie_out_partisanship_reason_ven.plot + 
  plot_layout(ncol=2, widths=c(1,1))

@

<<int_conj_skill, eval = TRUE, echo = FALSE, tidy=TRUE, fig.width = 8, fig.height = 3.5, out.width= "1\\linewidth", fig.align='center', warning=FALSE, message=FALSE, strip.white=TRUE, fig.cap="Interaction AMCE between partisanship and skill.">>=

#### Colombians

## Run interactive model
cjoint_inter_col <- CausalANOVA(chosen_profile ~ Partisanship + `Identity of Migrant` + 
                       `Skill Level` + `Probability of Employment` +
                       `Reason for Leaving` + Race + Gender,
       int2.formula = ~Partisanship:`Skill Level`,
       data = df_cjoint_clean_col, pair.id = df_cjoint_clean_col$pair_id, diff = TRUE,
       cluster = as.character(df_cjoint_clean_col$id), nway = 2, verbose = FALSE)

## Effect of partisanship conditional on skill
ceff_pol_skill_col <- ConditionalEffect(cjoint_inter_col, 
                                     treat.fac = "Partisanship", 
                                     cond.fac = "`Skill Level`", 
                                     base.ind = 1, verbose = FALSE)

ceff_pol_skill_col <- bind_rows(
  ceff_pol_skill_col$ConditionalEffects %>% .[[1]] %>%  
    as.data.frame() %>%
    mutate(`Skill Level` = "Low", Partisanship = rownames(.)),
  ceff_pol_skill_col$ConditionalEffects %>% .[[2]] %>%  
    as.data.frame() %>%
    mutate(`Skill Level` = "Medium", Partisanship = rownames(.)),
  ceff_pol_skill_col$ConditionalEffects %>% .[[3]] %>%  
    as.data.frame() %>%
    mutate(`Skill Level` = "High", 
           Partisanship = rownames(.)), 
) %>% filter(Partisanship != "Left")

ceff_pol_skill_col$Partisanship <- fct_relevel(ceff_pol_skill_col$Partisanship, "Right", "Center")
ceff_pol_skill_col$`Skill Level` <- fct_relevel(ceff_pol_skill_col$`Skill Level`, "Low", "Medium", "High")

amcie_out_partisanship_skill_col.plot <- ggplot(ceff_pol_skill_col, 
                                            aes(Partisanship, 
                                                ConditionalEffect, 
                                                group = `Skill Level`, 
                                                colour = `Skill Level`)) + 
  geom_point(position = position_dodge(.4)) + 
  geom_errorbar(aes(ymin = `2.5% CI`, ymax = `97.5% CI`), width = 0, 
                size = 1, position = position_dodge(.4)) + 
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  ylim(-.05, .2) +
  ggtitle("Colombian respondents") +
  ylab("Change in Pr(Migrant Preferred to Stay\nin the City and receive Benefits)") +
  scale_colour_manual(values=cbPalette) +
  theme_bw() +
  theme(legend.justification=c(1,0), 
        legend.position="none",
        legend.title = element_blank(),
        legend.key.size = unit(0.1, "cm")) +
  coord_flip()


#### Venezuelans
## Run interactive model
cjoint_inter_ven <- CausalANOVA(chosen_profile ~ Partisanship + `Identity of Migrant` + 
                       `Skill Level` + `Probability of Employment` +
                       `Reason for Leaving` + Race + Gender,
       int2.formula = ~Partisanship:`Skill Level`,
       data = df_cjoint_clean_ven, pair.id = df_cjoint_clean_ven$pair_id, diff = TRUE,
       cluster = as.character(df_cjoint_clean_ven$id), nway = 2, verbose = FALSE)

## Effect of partisanship conditional on skill
ceff_pol_skill_ven <- ConditionalEffect(cjoint_inter_ven, 
                                     treat.fac = "Partisanship", 
                                     cond.fac = "`Skill Level`", 
                                     base.ind = 1, verbose = FALSE)

ceff_pol_skill_ven <- bind_rows(
  ceff_pol_skill_ven$ConditionalEffects %>% .[[1]] %>%  
    as.data.frame() %>%
    mutate(`Skill Level` = "Low", Partisanship = rownames(.)),
  ceff_pol_skill_ven$ConditionalEffects %>% .[[2]] %>%  
    as.data.frame() %>%
    mutate(`Skill Level` = "Medium", Partisanship = rownames(.)),
  ceff_pol_skill_ven$ConditionalEffects %>% .[[3]] %>%  
    as.data.frame() %>%
    mutate(`Skill Level` = "High", 
           Partisanship = rownames(.)), 
) %>% filter(Partisanship != "Left")

ceff_pol_skill_ven$Partisanship <- fct_relevel(ceff_pol_skill_col$Partisanship, "Right", "Center")
ceff_pol_skill_ven$`Skill Level` <- fct_relevel(ceff_pol_skill_col$`Skill Level`, "Low", "Medium", "High")

amcie_out_partisanship_skill_ven.plot <- ggplot(ceff_pol_skill_ven, 
                                            aes(Partisanship, 
                                                ConditionalEffect, 
                                                group = `Skill Level`, 
                                                colour = `Skill Level`)) + 
  geom_point(position = position_dodge(.4)) + 
  geom_errorbar(aes(ymin = `2.5% CI`, ymax = `97.5% CI`), width = 0, 
                size = 1, position = position_dodge(.4)) + 
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  ylim(-.05, .2) +
  ggtitle("Venezuelan respondents") +
  ylab("Change in Pr(Migrant Preferred to Stay\nin the City and receive Benefits)") +
  scale_colour_manual(values=cbPalette) +
  theme_bw() +
  theme(legend.justification=c(1,0), 
        legend.position=c(.99,0.01),
        legend.title = element_blank(),
        legend.key.size = unit(0.1, "cm")) +
  coord_flip()

## plot
amcie_out_partisanship_skill_col.plot + 
  amcie_out_partisanship_skill_ven.plot + 
  plot_layout(ncol=2, widths=c(1,1))

@

<<int_conj_employability, eval = TRUE, echo = FALSE, tidy=TRUE, fig.width = 8, fig.height = 3.5, out.width= "1\\linewidth", fig.align='center', warning=FALSE, message=FALSE, strip.white=TRUE, fig.cap="Interaction AMCE between partisanship and employability">>=

#### Colombians

## Run interactive model
cjoint_inter_col <- CausalANOVA(chosen_profile ~ Partisanship + `Identity of Migrant` + 
                       `Skill Level` + `Probability of Employment` +
                       `Reason for Leaving` + Race + Gender,
       int2.formula = ~Partisanship:`Probability of Employment`,
       data = df_cjoint_clean_col, pair.id = df_cjoint_clean_col$pair_id, diff = TRUE,
       cluster = as.character(df_cjoint_clean_col$id), nway = 2, verbose = FALSE)

## Effect of partisanship conditional on employment
ceff_pol_employ_col <- ConditionalEffect(cjoint_inter_col, 
                                     treat.fac = "Partisanship", 
                                     cond.fac = "`Probability of Employment`", 
                                     base.ind = 1, verbose = FALSE)

ceff_pol_employ_col <- bind_rows(
  ceff_pol_employ_col$ConditionalEffects %>% .[[1]] %>%  
    as.data.frame() %>%
    mutate(`Probability of Employment` = "Cannot work", Partisanship = rownames(.)),
  ceff_pol_employ_col$ConditionalEffects %>% .[[2]] %>%  
    as.data.frame() %>%
    mutate(`Probability of Employment` = "Unlikely to find work", Partisanship = rownames(.)),
  ceff_pol_employ_col$ConditionalEffects %>% .[[3]] %>%  
    as.data.frame() %>%
    mutate(`Probability of Employment` = "Very likely to find work", 
           Partisanship = rownames(.)), 
) %>% filter(Partisanship != "Left")

ceff_pol_employ_col$Partisanship <- fct_relevel(ceff_pol_employ_col$Partisanship, "Right", "Center")
ceff_pol_employ_col$`Probability of Employment` <- fct_relevel(ceff_pol_employ_col$`Probability of Employment`, 
                                                              "Cannot work", "Unlikely to find work", 
                                                              "Very likely to find work")

amcie_out_partisanship_employ_col.plot <- ggplot(ceff_pol_employ_col, 
                                            aes(Partisanship, 
                                                ConditionalEffect, 
                                                group = `Probability of Employment`, 
                                                colour = `Probability of Employment`)) + 
  geom_point(position = position_dodge(.4)) + 
  geom_errorbar(aes(ymin = `2.5% CI`, ymax = `97.5% CI`), width = 0, 
                size = 1, position = position_dodge(.4)) + 
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  ylim(-.05, .2) +
  ggtitle("Colombian respondents") +
  ylab("Change in Pr(Migrant Preferred to Stay\nin the City and receive Benefits)") +
  scale_colour_manual(values=cbPalette) +
  theme_bw() +
  theme(legend.justification=c(1,0), 
        legend.position="none",
        legend.title = element_blank(),
        legend.key.size = unit(0.1, "cm")) +
  coord_flip()


#### Venezuelans
## Run interactive model
cjoint_inter_ven <- CausalANOVA(chosen_profile ~ Partisanship + `Identity of Migrant` + 
                       `Skill Level` + `Probability of Employment` +
                       `Reason for Leaving` + Race + Gender,
       int2.formula = ~Partisanship:`Probability of Employment`,
       data = df_cjoint_clean_ven, pair.id = df_cjoint_clean_ven$pair_id, diff = TRUE,
       cluster = as.character(df_cjoint_clean_ven$id), nway = 2, verbose = FALSE)

## Effect of partisanship conditional on reason for leaving
ceff_pol_employ_ven <- ConditionalEffect(cjoint_inter_ven, 
                                     treat.fac = "Partisanship", 
                                     cond.fac = "`Probability of Employment`", 
                                     base.ind = 1, verbose = FALSE)

ceff_pol_employ_ven <- bind_rows(
  ceff_pol_employ_ven$ConditionalEffects %>% .[[1]] %>%  
    as.data.frame() %>%
    mutate(`Probability of Employment` = "Cannot work", Partisanship = rownames(.)),
  ceff_pol_employ_ven$ConditionalEffects %>% .[[2]] %>%  
    as.data.frame() %>%
    mutate(`Probability of Employment` = "Unlikely to find work", Partisanship = rownames(.)),
  ceff_pol_employ_ven$ConditionalEffects %>% .[[3]] %>%  
    as.data.frame() %>%
    mutate(`Probability of Employment` = "Very likely to find work", 
           Partisanship = rownames(.)), 
) %>% filter(Partisanship != "Left")

ceff_pol_employ_ven$Partisanship <- fct_relevel(ceff_pol_employ_ven$Partisanship, "Right", "Center")
ceff_pol_employ_ven$`Probability of Employment` <- fct_relevel(ceff_pol_employ_ven$`Probability of Employment`, 
                                                              "Cannot work", "Unlikely to find work", 
                                                              "Very likely to find work")


amcie_out_partisanship_employ_ven.plot <- ggplot(ceff_pol_employ_ven, 
                                            aes(Partisanship, 
                                                ConditionalEffect, 
                                                group = `Probability of Employment`, 
                                                colour = `Probability of Employment`)) + 
  geom_point(position = position_dodge(.4)) + 
  geom_errorbar(aes(ymin = `2.5% CI`, ymax = `97.5% CI`), width = 0, 
                size = 1, position = position_dodge(.4)) + 
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  ylim(-.05, .2) +
  ggtitle("Venezuelan respondents") +
  ylab("Change in Pr(Migrant Preferred to Stay\nin the City and receive Benefits)") +
  scale_colour_manual(values=cbPalette) +
  theme_bw() +
  theme(legend.justification=c(1,0), 
        legend.position=c(.99,0.01),
        legend.title = element_blank(),
        legend.key.size = unit(0.1, "cm")) +
  coord_flip()

## plot
amcie_out_partisanship_employ_col.plot + 
  amcie_out_partisanship_employ_ven.plot + 
  plot_layout(ncol=2, widths=c(1,1))


@

\newpage
\section{Marginal Means for Conditional Conjoint Results}
\label{SIsec:mm}

<<conj_prep, eval = TRUE, echo = FALSE, tidy=TRUE>>=

#### Colombians

## Full conjoint dataframe
df_cjoint_col <- df_col %>% 
  dplyr::select(matches("con[1-5]"), 
                starts_with("rd_"),
                partisanship,
                skilled_labor,
                salary,
                benefits_bi,
                indir_contact_index_2,
                city,
                ven_friends,
                ven_beg_bi,
                race) %>% 
  dplyr::mutate(id = 1:n(),
                partisanship = as.factor(partisanship),
                skilled_labor = as.factor(skilled_labor),
                salary = as.factor(salary),
                benefits_bi = as.factor(benefits_bi),
                indir_contact_index_2 = as.factor(indir_contact_index_2),
                city = as.factor(city),
                ven_friends = as.factor(ven_friends),
                ven_beg_bi = as.factor(ven_beg_bi),
                race = as.factor(race))

## Clean responses
df_cjoint_response_col <- df_cjoint_col %>%
  dplyr::select(id, matches("con[1-5]")) %>%
  gather(key = round, value = profile_choice, -id) %>%
  mutate(round = parse_number(round))

## Clean attributes
df_cjoint_covs_col <- df_cjoint_col %>%
  dplyr::select(id, starts_with("rd_")) %>%
  gather(key = feature_round_profile, value = attribute, -id) %>%
  mutate(round = parse_number(feature_round_profile),
         profile = case_when(grepl("_a_", feature_round_profile)~"a",
                             grepl("_b_", feature_round_profile)~"b",
                             TRUE~"ruh-roh"),
         feature = gsub("rd_[1-5]_[a|b]_", "", feature_round_profile)) %>%
  dplyr::select(-feature_round_profile) %>%
  spread(key = feature, value = attribute)

## Bind responses and attributes together
df_cjoint_clean_col <- inner_join(df_cjoint_response_col, df_cjoint_covs_col,
                              by = c("id", "round")) %>%
  left_join(df_cjoint_col %>% dplyr::select(id, 
                                            partisanship,
                                            skilled_labor,
                                            salary,
                                            benefits_bi,
                                            indir_contact_index_2,
                                            city,
                                            ven_friends,
                                            ven_beg_bi,
                                            race)) %>%
  mutate(chosen_profile = case_when(profile_choice == 1 & profile == "a"~1,
                                    profile_choice == 2 & profile == "b"~1,
                                    profile_choice == 1 & profile == "b"~0,
                                    profile_choice == 2 & profile == "a"~0,
                                    TRUE~NA_real_),
         Partisanship = as.factor(Partisanship), 
         `Skill Level` = as.factor(`Skill Level`), 
         `Identity of Migrant` = as.factor(`Identity of Migrant`),
         `Probability of Employment` = as.factor(`Probability of Employment`), 
         Race = as.factor(Race), 
         `Reason for Leaving` = as.factor(`Reason for Leaving`), 
         Gender = as.factor(Gender)
         ) %>%
  drop_na(-c(partisanship,
                skilled_labor,
                salary,
                benefits_bi,
                indir_contact_index_2,
                city,
                ven_friends,
                ven_beg_bi,
                race)) %>%
  mutate(pair_id = paste0(id, round)) %>%
  group_by(pair_id) %>%
  filter(n() == 2) %>% ungroup()

## Run conjoint estimator
set_baseline <- list(Partisanship = "Left",
                     `Skill Level` = "Low",
                     `Identity of Migrant` = "Colombians displaced from around Cali",
                     `Probability of Employment` = "Cannot work",
                      Race = "Mestizo",
                     `Reason for Leaving` = "Fear of crime in their area")

df_cjoint_clean_col$Partisanship <- fct_relevel(df_cjoint_clean_col$Partisanship, "Left", "Center", "Right")
df_cjoint_clean_col$`Skill Level` <- fct_relevel(df_cjoint_clean_col$`Skill Level`, "Low", "Medium", "High")

#### Venezuelans

## Full conjoint dataframe
df_cjoint_ven <- df_ven %>% 
 dplyr::select(matches("con[1-5]"), 
        starts_with("rd_"),
        partisanship,
        skilled_labor,
        salary,
        benefits_bi,
        city) %>% 
 dplyr::mutate(id = 1:n(),
        partisanship = as.factor(partisanship),
        skilled_labor = as.factor(skilled_labor),
        salary = as.factor(salary),
        benefits_bi = as.factor(benefits_bi),
        city = as.factor(city))

## Clean responses
df_cjoint_response_ven <- df_cjoint_ven %>%
 dplyr::select(id, matches("con[1-5]")) %>%
 gather(key = round, value = profile_choice, -id) %>%
 mutate(round = parse_number(round))

## Clean attributes
df_cjoint_covs_ven <- df_cjoint_ven %>%
 dplyr::select(id, starts_with("rd_")) %>%
 gather(key = feature_round_profile, value = attribute, -id) %>%
 mutate(round = parse_number(feature_round_profile),
     profile = case_when(grepl("_a_", feature_round_profile)~"a",
               grepl("_b_", feature_round_profile)~"b",
               TRUE~"ruh-roh"),
     feature = gsub("rd_[1-5]_[a|b]_", "", feature_round_profile)) %>%
 dplyr::select(-feature_round_profile) %>%
 spread(key = feature, value = attribute)

## Bind responses and attributes together
df_cjoint_clean_ven <- inner_join(df_cjoint_response_ven, df_cjoint_covs_ven,
               by = c("id", "round")) %>%
 left_join(df_cjoint_ven %>% dplyr::select(id,
                      partisanship,
                      skilled_labor,
                      salary,
                      benefits_bi,
                      city)) %>%
 mutate(chosen_profile = case_when(profile_choice == 1 & profile == "a"~1,
                  profile_choice == 2 & profile == "b"~1,
                  profile_choice == 1 & profile == "b"~0,
                  profile_choice == 2 & profile == "a"~0,
                  TRUE~NA_real_),
     Partisanship = as.factor(Partisanship), 
     `Skill Level` = as.factor(`Skill Level`), 
     `Identity of Migrant` = as.factor(`Identity of Migrant`),
     `Probability of Employment` = as.factor(`Probability of Employment`), 
     Race = as.factor(Race), 
     `Reason for Leaving` = as.factor(`Reason for Leaving`), 
     Gender = as.factor(Gender)
     ) %>%
 drop_na(-c(partisanship,
        skilled_labor,
        salary,
        benefits_bi,
        city)) %>%
 mutate(pair_id = paste0(id, round)) %>%
 group_by(pair_id) %>%
 filter(n() == 2) %>% ungroup()
 #drop_na(profile_choice)

## Run conjoint estimator
set_baseline <- list(Partisanship = "Left",
           `Skill Level` = "Low",
           `Identity of Migrant` = "Colombians displaced from around Cali",
           `Probability of Employment` = "Cannot work",
           Race = "Mestizo",
           `Reason for Leaving` = "Fear of crime in their area;")

df_cjoint_clean_ven$Partisanship <- fct_relevel(df_cjoint_clean_ven$Partisanship, "Left", "Center", "Right")
df_cjoint_clean_ven$`Skill Level` <- fct_relevel(df_cjoint_clean_ven$`Skill Level`, "Low", "Medium", "High")


@

This section shows the marginal means estimates for Figure \ref{fig:cond_conj} in the paper using the procedure described in \citet{leeper2020measuring}. Whereas the ACME conditional on respondent polarization and city, respectively, in those figures show the causal effect of the features shown relative to a baseline attribute, comparison across the subgroups (e.g. respondents who are Left, Center, Right) may be misleading if those subgroups diverge in preferences towards the baseline category. Thus, the marginal mean shows the average for every level of the attribute of interest, marginalizing over all over attributes. This estimate helps us descriptively compare across subgroups, and allows us to check whether subgroups do differ on there preference for the baseline category. Figures \ref{fig:cond_conj_part_mm} and \ref{fig:cond_conj_city_mm} below confirm that our interpretations comparing across subgroups using AMCEs hold.

<<cond_conj_part_mm, eval = TRUE, echo = FALSE, tidy=TRUE, fig.width = 7, fig.height = 6, out.width= ".8\\linewidth", fig.align='center', warning=FALSE, message=FALSE, strip.white=TRUE, fig.cap="Marginal means of migrant profile political polarization conditional on respondent political polarization.">>=

#### Colombians
col_mm_uncond <- cregg::cj(
  chosen_profile ~ Partisanship + `Identity of Migrant` + `Skill Level` +
    `Probability of Employment` +`Reason for Leaving` + Race + Gender,
  data = df_cjoint_clean_col, 
  estimate = "mm",
  id = ~id
)
col_mm_condpart <- cregg::cj(
  chosen_profile ~ Partisanship + `Identity of Migrant` + `Skill Level` +
    `Probability of Employment` +`Reason for Leaving` + Race + Gender,
  data = df_cjoint_clean_col, 
  estimate = "mm",
  id = ~id,
  by = ~partisanship
)

#### Venezuelans
ven_mm_uncond <- cregg::cj(
  chosen_profile ~ Partisanship + `Identity of Migrant` + `Skill Level` +
    `Probability of Employment` +`Reason for Leaving` + Race + Gender,
  data = df_cjoint_clean_ven[!is.na(df_cjoint_clean_ven$partisanship),], 
  estimate = "mm",
  id = ~id
)
ven_mm_condpart <- cregg::cj(
  chosen_profile ~ Partisanship + `Identity of Migrant` + `Skill Level` +
    `Probability of Employment` +`Reason for Leaving` + Race + Gender,
  data = df_cjoint_clean_ven[!is.na(df_cjoint_clean_ven$partisanship),], 
  estimate = "mm",
  id = ~id,
  by = ~partisanship
)

## Create plotting dataframe
plot_col <- col_mm_uncond %>%
  filter(feature == "Partisanship") %>%
  mutate(partisanship = "Unconditional") %>%
  bind_rows(
    col_mm_condpart %>%
      filter(feature == "Partisanship") %>%
      mutate(partisanship = as.character(partisanship),
             partisanship = case_when(
               partisanship == "1"~"Conditional on Left",
               partisanship == "2"~"Conditional on Center",
               partisanship == "3"~"Conditional on Right",
               TRUE~"ruh-roh"
              ))
  ) %>%
  mutate(partisanship = fct_relevel(partisanship,
                                    "Unconditional",
                                    "Conditional on Left",
                                    "Conditional on Center",
                                    "Conditional on Right"
                                    ),
         level = fct_relevel(level, "Right", "Center", "Left"))

plot_ven <- ven_mm_uncond %>%
  filter(feature == "Partisanship") %>%
  mutate(partisanship = "Unconditional") %>%
  bind_rows(
    ven_mm_condpart %>%
      filter(feature == "Partisanship") %>%
      mutate(partisanship = as.character(partisanship),
             partisanship = case_when(
               partisanship == "1"~"Conditional on Left",
               partisanship == "2"~"Conditional on Center",
               partisanship == "3"~"Conditional on Right",
               TRUE~"ruh-roh"
              ))
  ) %>%
  mutate(partisanship = fct_relevel(partisanship,
                                    "Unconditional",
                                    "Conditional on Left",
                                    "Conditional on Center",
                                    "Conditional on Right"
                                    ),
         level = fct_relevel(level, "Right", "Center", "Left"))

## Plot
mm_col_plot <- ggplot(plot_col, aes(level, estimate)) + 
  geom_point(size = 2.5, colour = "#D55E00") + 
  geom_errorbar(aes(ymin = lower, ymax = upper), 
                width = 0, size = .5, colour = "#D55E00") + 
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  #ylim(-.15, .2) +
  facet_wrap(~partisanship, ncol = 1) + 
  ggtitle("Colombian respondents") +
  ylab("Marginal Mean of Pr(Preferred Migrant)") +
  theme_bw() + 
  theme(legend.position = "none",
        axis.title.y = element_blank()) +
  coord_flip()

mm_ven_plot <- ggplot(plot_ven, aes(level, estimate)) + 
  geom_point(size = 2.5, colour = "#D55E00") + 
  geom_errorbar(aes(ymin = lower, ymax = upper), 
                width = 0, size = .5, colour = "#D55E00") + 
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  #ylim(-.15, .2) +
  facet_wrap(~partisanship, ncol = 1) + 
  ggtitle("Venezuelan respondents") +
  ylab("Marginal Mean of Pr(Preferred Migrant)") +
  theme_bw() + 
  theme(legend.position = "none",
        axis.title.y = element_blank()) +
  coord_flip()

#plot 
mm_col_plot + mm_ven_plot + plot_layout(ncol=2, widths=c(1,1))

@

<<cond_conj_city_mm, eval = TRUE, echo = FALSE, tidy=TRUE, fig.width = 7, fig.height = 6, out.width= ".8\\linewidth", fig.align='center', warning=FALSE, message=FALSE, strip.white=TRUE, fig.cap="Marginal means of migrant profile political polarization conditional on respondent's survey city.">>=

#### Colombians
col_mm_uncond <- cregg::cj(
  chosen_profile ~ Partisanship + `Identity of Migrant` + `Skill Level` +
    `Probability of Employment` +`Reason for Leaving` + Race + Gender,
  data = df_cjoint_clean_col, 
  estimate = "mm",
  id = ~id
)
col_mm_condpart <- cregg::cj(
  chosen_profile ~ Partisanship + `Identity of Migrant` + `Skill Level` +
    `Probability of Employment` +`Reason for Leaving` + Race + Gender,
  data = df_cjoint_clean_col, 
  estimate = "mm",
  id = ~id,
  by = ~city
)

#### Venezuelans
ven_mm_uncond <- cregg::cj(
  chosen_profile ~ Partisanship + `Identity of Migrant` + `Skill Level` +
    `Probability of Employment` +`Reason for Leaving` + Race + Gender,
  data = df_cjoint_clean_ven[!is.na(df_cjoint_clean_ven$city),], 
  estimate = "mm",
  id = ~id
)
ven_mm_condpart <- cregg::cj(
  chosen_profile ~ Partisanship + `Identity of Migrant` + `Skill Level` +
    `Probability of Employment` +`Reason for Leaving` + Race + Gender,
  data = df_cjoint_clean_ven[!is.na(df_cjoint_clean_ven$city),], 
  estimate = "mm",
  id = ~id,
  by = ~city
)

## Create plotting dataframe
plot_col <- col_mm_uncond %>%
  filter(feature == "Partisanship") %>%
  mutate(city = "Unconditional") %>%
  bind_rows(
    col_mm_condpart %>%
      filter(feature == "Partisanship") %>%
      mutate(city = as.character(city),
             city = case_when(
               city == "Cali"~"Conditional on Cali",
               city == "Cúcuta"~"Conditional on Cúcuta",
               TRUE~"ruh-roh"
              ))
  ) %>%
     mutate(city = fct_relevel(city, "Unconditional"),
            level = fct_relevel(level, "Right", "Center", "Left"))

plot_ven <- ven_mm_uncond %>%
  filter(feature == "Partisanship") %>%
  mutate(city = "Unconditional") %>%
  bind_rows(
    ven_mm_condpart %>%
      filter(feature == "Partisanship") %>%
      mutate(city = as.character(city),
             city = case_when(
               city == "Cali"~"Conditional on Cali",
               city == "Cúcuta"~"Conditional on Cúcuta",
               TRUE~"ruh-roh"
              ))
  ) %>%
  mutate(city = fct_relevel(city, "Unconditional"),
         level = fct_relevel(level, "Right", "Center", "Left"))

## Plot
mm_col_plot <- ggplot(plot_col, aes(level, estimate)) + 
  geom_point(size = 2.5, colour = "#D55E00") + 
  geom_errorbar(aes(ymin = lower, ymax = upper), 
                width = 0, size = .5, colour = "#D55E00") +  
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  #ylim(-.15, .2) +
  facet_wrap(~city, ncol = 1) + 
  ggtitle("Colombian respondents") +
  ylab("Marginal Mean of Pr(Preferred Migrant)") +
  theme_bw() + 
  theme(legend.position = "none",
        axis.title.y = element_blank()) +
  coord_flip()

mm_ven_plot <- ggplot(plot_ven, aes(level, estimate)) + 
  geom_point(size = 2.5, colour = "#D55E00") + 
  geom_errorbar(aes(ymin = lower, ymax = upper), 
                width = 0, size = .5, colour = "#D55E00") + 
  geom_hline(aes(yintercept = 0), lty = "dashed") + 
  #ylim(-.15, .2) +
  facet_wrap(~city, ncol = 1) + 
  ggtitle("Venezuelan respondents") +
  ylab("Marginal Mean of Pr(Preferred Migrant)") +
  theme_bw() + 
  theme(legend.position = "none",
        axis.title.y = element_blank()) +
  coord_flip()

#plot 
mm_col_plot + mm_ven_plot + plot_layout(ncol=2, widths=c(1,1))

@


\newpage
\section{Tweets by \'Alvaro Uribe mobilizing voters for the 2018 Elections}
\label{SIsec:tweets}

This section presents examples of former President \'Alvaro Uribe's tweets in the lead-up to Colombia's presidential elections on 27 May 2018. One of the most common themes was a concern that Colombia would become a ``second Venezuela'' followed by a call to vote for Iv\'an Duque:

\begin{figure}[H]
\centering
\begin{minipage}{.3\textwidth}
  \centering
  \includegraphics[width=1\linewidth]{tweet.Dec3}
  \captionof{figure}{\scriptsize{@AlvaroUribeVel. ``No vamos a dejar que a este país lo vuelvan una segunda Venezuela. @CeDemocratico la esperanza. Nos mueve el amor por Colombia. Con Mano firme y corazón grande, retomaremos el rumbo. Nuestra obsesión: cero corrupción, bajar impuestos, subir salarios'' Twitter, 3 Dec. 2017, \url{twitter.com/AlvaroUribeVel/status/937480548271972353?s=20}.}}
\end{minipage}% 
\hspace{.4cm}%
\begin{minipage}{.3\textwidth}
  \centering
  \includegraphics[width=1\linewidth]{tweet.Jan25}
  \captionof{figure}{\scriptsize{@AlvaroUribeVel. ``En las calles de Armenia, el Quindío, pedacito de Cielo, no podemos ser una segunda Venezuela, con Iván Duque y nuestra coalición por una democracia con desempeño sobresaliente en seguridad, emprendimiento y equidad.'' Twitter, 25 Jan 2018, \url{twitter.com/AlvaroUribeVel/status/956636838768398339?s=20}.}}
\end{minipage}%
\hspace{.4cm}%
\begin{minipage}{.3\textwidth}
  \centering
  \includegraphics[width=1\linewidth]{tweet.Feb8}
  \captionof{figure}{\scriptsize{@AlvaroUribeVel. ``Por un país seguro, con empresas vigorosas, con menos impuestos, con trabajadores reivindicados, con mejores salarios; un país solidario. \mbox{!`}No a la lucha de clases! \mbox{!`}No a una segunda Venezuela!'' Twitter, 8 Feb. 2018, \url{twitter.com/AlvaroUribeVel/status/961572834685276160?s=20}.}}
\end{minipage}
\end{figure}

\begin{figure}[H]
\centering
\begin{minipage}{.3\textwidth}
  \centering
  \includegraphics[width=1\linewidth]{tweet.Feb21}
  \captionof{figure}{\scriptsize{@AlvaroUribeVel. ``Votemos por Colombia: -Estado austero, sin corrupción -Empresas vigorosas con menos impuestos -Trabajadores reivindicados con mejores salarios - 6 días al año sin IVA para productos básicos -Pais solidario sin odio de clases -No a una 2a Venezuela Iván Duque Presidente.'' Twitter, 21 Feb 2018, \url{twitter.com/AlvaroUribeVel/status/966478804347686912?s=20}.}}
\end{minipage}%
\hspace{.4cm}%
\begin{minipage}{.3\textwidth}
  \centering
  \includegraphics[width=1\linewidth]{tweet.Feb24}
  \captionof{figure}{\scriptsize{@AlvaroUribeVel. ``Ciudadanos, escuchen los audios a las interceptaciones a mis llamadas, no he delinquido, he procurado desmontar la violencia moral de los narco terroristas, sus voceros y aliados. No a una segunda Venezuela. Ayúdennos a ganar.'' Twitter, 24 Feb 2018, \url{twitter.com/AlvaroUribeVel/status/967599110428733442?s=20}.}}
\end{minipage}%
\hspace{.4cm}%
\begin{minipage}{.3\textwidth}
  \centering
  \includegraphics[width=1\linewidth]{tweet.Feb26}
  \captionof{figure}{\scriptsize{@AlvaroUribeVel. ``Gracias Cartagena,pueblo heroico sabe decir no al odio de clases. No segunda Venezuela Cero corrupción,cero dosis personal. Patria segura con amor a los soldados y policías.Sin impunidad. Menos impuestos más salario, 6 días x año sin IVA a lo básico. País cristiano,solidario.'' Twitter, 26 Feb. 2018, \url{twitter.com/AlvaroUribeVel/status/968267406618644480?s=20}.}}
\end{minipage}
\end{figure}

\begin{figure}[H]
\centering
\begin{minipage}{.35\textwidth}
  \centering
  \includegraphics[width=1\linewidth]{tweet.Apr14}
  \captionof{figure}{\scriptsize{@AlvaroUribeVel. ``Funza, el deterioro de Colombia no puede servir de excusa para llevarnos a ser una segunda Venezuela. Iv\'an Duque Presidente. Marta Lucía Ramírez VicePte.'' Twitter, 14 Apr. 2018, \url{twitter.com/AlvaroUribeVel/status/985288392161325057?s=20}.}}
\end{minipage}% 
\hspace{.4cm}%
\begin{minipage}{.35\textwidth}
  \centering
  \includegraphics[width=1\linewidth]{tweet.May6}
  \captionof{figure}{\scriptsize{@AlvaroUribeVel. ``Gracias Montería, el odio de clases empeora más el deteriora y lleva a una segunda Venezuela. Por un país fraterno. Iván Duque Pte Marta. Lucía Ramírez, VicePte.'' Twitter, 6 May 2018, \url{twitter.com/AlvaroUribeVel/status/993266357398536197?s=20}.}}
\end{minipage}
\end{figure}

Uribe's tweets also focused on threats from socialism or Castro-Chavismo coming to Colombia. For example:

\begin{figure}[H]
\centering
\begin{minipage}{.35\textwidth}
  \centering
  \includegraphics[width=1\linewidth]{tweet.Feb6}
  \captionof{figure}{\scriptsize{@AlvaroUribeVel. ``Gracias Floridablanca: “no permitamos un socialismo” me dicen compatriotas de Santander y hermanos de Venezuela'' Twitter, 6 Feb. 2018, \url{twitter.com/AlvaroUribeVel/status/963524737992413185?s=20}.}}
\end{minipage}% 
\hspace{.4cm}%
\begin{minipage}{.35\textwidth}
  \centering
  \includegraphics[width=1\linewidth]{tweet.Feb13}
  \captionof{figure}{\scriptsize{@AlvaroUribeVel. ``Ni los colombianos ni la comunidad internacional alcanzamos a imaginar el poder criminal de Farc y Eln (cobalto) en Venezuela, en apoyo a la tiranía. La paz ha sido la legalización de algunos para que los otros sostengan al Castro Chavismo en Venezuela y lo impongan en Colombia'' Twitter, 13 Feb. 2018, \url{twitter.com/AlvaroUribeVel/status/960851461046448129?s=20}.}}
\end{minipage}
\end{figure}


\newpage
\section{Observational Results: Factors Associated with Openness Toward Migrants}
\label{SIsec:regress}

In this section, we present the results with observational data for each of the hypotheses that we outlined in our pre-analysis plan and reproduce below. To be clear, we are not making causal claims with this observational data; but we are presenting descriptive associations that are important for understanding what types of Colombians tend to feel more inclusive or exclusionary towards Venezuelan migrants.

We include a set of descriptive questions to understand Colombians' broader anxieties about Venezuelan migration. These questions examine the major concerns identified in research on immigrant reception. Respondents were given several questions on which they could agree or disagree.\footnote{Within each block, the statements were randomized to either be positive or negative to reduce wording effects. We find no evidence of wording effects.} Our main observational outcome of interest is {\em{Openness}}, which is an index to capture the extent to which Colombians approve of hosting Venezuelans and allowing additional migrants to enter. We sum over the following questions based on the scores shown, and then we rescale the index to vary between 0 and 1 for interpretability.

\begin{singlespace}
\vspace{-.5cm}
\texttt{
\begin{enumerate}
%\setlength{\itemsep}{-10pt}
  \item There are too many Venezuelans in Colombia. Agree = 0, Disagree = 1
  \item Colombia should close its borders immediately; it isn't possible to accept more migrants at this moment. Agree = 0, Disagree = 1
  \item Which of the following statements aligns closest with what you believe Colombia should do about Venezuelan migrants:
  \begin{itemize}
   %\setlength{\itemsep}{-10pt}
   \item Colombia should send back the majority of Venezuelans to Venezuela now = 0
   \item Colombia should send back the majority of Venezuelans to Venezuela once there has been a political transition =.5
   \item Colombia should send the majority of Venezuelans to Venezuela once the country has recovered economically = .5
   \item Colombia should let Venezuelans who desire to stay in Colombia permanently = 1
\end{itemize}
\end{enumerate}
}
\end{singlespace}


We expect {\em{O1. Colombians who hold political misperceptions are less likely to support open border policies.}}.  It is possible that respondents hold political misperceptions because they want to close the border for unrelated reasons. We therefore also look at the broader determinants of attitudes towards open migration politics to test alternative explanations. We regress our key dependent variable on the relevant independent variables using OLS. In these models, we control for respondent demographics: city, gender, age, race, education level, number of children, marital status, religion, religiosity, wealth (based on ownership of 14 household items), skill level, salary level (4-point scale), political views, and access to welfare benefits (based on the household's receipt of the five most common social benefits). We expect political misperceptions to be more important than other common concerns in predicting border attitudes. 

For ease of interpreting coefficients, we rescale the Openness Index outcome variable to vary between minimum 0 and maximum 1, as well as the independent variables measuring: left-center-right partisanship, skilled labor index, contract formality index, salary measure, benefits index, direct contact with Venezuelans index, indirect contact with Venezuelans index, and cultural similarity index.

<<mainreg, eval = TRUE, echo = FALSE, tidy=TRUE, fig.width = 5, fig.height = 5.5, out.width= ".8\\linewidth", fig.align='center', warning=FALSE, message=FALSE, strip.white=TRUE, fig.cap="OLS estimates with 95\\% CIs regressing Openness to Venezuelans (Index) onto variables mapping onto political ideology (red), labor market (blue), fiscal (green), sociotropic (black), and cultural (yellow) concerns. Openness Index and all independent variables shown on the left have been rescaled to vary between 0 and 1.">>=

### Main regression analysis fig for the paper

# Only Openness index as outcomes, removing some of the attitude on attitude assocs
open_main_lm <- df_col %>% 
 dplyr::select(
  open_index_res, 
  pol_left, 
  pol_guerilla,
  partisanship_res,
  skilled_labor_res, 
  contract_res, 
  salary_res, 
  benefits_index_res, 
  dir_contact_index_res, 
  indir_contact_index_res, 
  int_disp, 
  int_disp_fam, 
  city, 
  hist_culture,
  male, 
  age, 
  race, 
  education, 
  kids, 
  marriage, 
  religion2, 
  religiosity, 
  wealth_index
 ) %>%
 pivot_longer(cols = c(pol_left, pol_guerilla, contract_res, 
            dir_contact_index_res, indir_contact_index_res, 
            int_disp, int_disp_fam, hist_culture),
        names_to = "indep_var") %>%
 group_by(indep_var) %>%
 do({
  lm_robust(update(open_index_res ~ value, reformulate(c(".", labels(demo_ctrls_col)))), data = .) %>% 
   tidy() %>%
   filter(term == "value")
 }) %>% bind_rows(
  lm_robust(update(open_index_res ~ 1, reformulate(c(".", labels(demo_ctrls_col)))), data = df_col) %>%
   tidy() %>%
   filter(term %in% c("cityCúcuta", "partisanship_res", "salary_res", "benefits_index_res", "skilled_labor_res")) %>%
   rename(indep_var = term) %>%
   mutate(term = "value")
 )

cbPalette2 <- c("#D55E00", "#0072B2", "#009E73", "#000000", "#E69F00")

open_main_lm.plot <- open_main_lm %>%
    as.data.frame() %>%
    mutate(order = case_when(indep_var == "partisanship_res"~13, #reordering for the plot
            indep_var == "pol_left"~12,
            indep_var == "pol_guerilla"~11,
            indep_var == "skilled_labor_res"~10,
            indep_var == "contract_res"~9,
            indep_var == "salary_res"~8,
            indep_var == "benefits_index_res"~7,
            indep_var == "cityCúcuta"~6,
            indep_var == "dir_contact_index_res"~5,
            indep_var == "indir_contact_index_res"~4,
            indep_var == "int_disp"~3,
            indep_var == "int_disp_fam"~2,
            indep_var == "hist_culture"~1,
            TRUE ~ NA_real_),
        group = case_when(indep_var == "partisanship_res"~"a_partisanship", #subgroup labels
            indep_var == "pol_left"~"a_partisanship",
            indep_var == "pol_guerilla"~"a_partisanship",
            indep_var == "skilled_labor_res"~"b_labor",
            indep_var == "contract_res"~"b_labor",
            indep_var == "salary_res"~"c_fiscal",
            indep_var == "benefits_index_res"~"c_fiscal",
            indep_var == "cityCúcuta"~"d_socio",
            indep_var == "dir_contact_index_res"~"e_cultural",
            indep_var == "indir_contact_index_res"~"e_cultural",
            indep_var == "int_disp"~"e_cultural",
            indep_var == "int_disp_fam"~"e_cultural",
            indep_var == "hist_culture"~"e_cultural",
            TRUE ~ NA_character_)) %>%
 ggplot(aes(x = fct_reorder(indep_var, order), 
       y = estimate, 
       colour = group, 
       group = group)) + 
 geom_point(size = 2.5) + 
 geom_errorbar(aes(ymin = conf.low, ymax = conf.high), 
        width = 0, size = .5) + 
 scale_colour_manual(values = cbPalette2) + 
 geom_hline(aes(yintercept = 0), lty = "dashed") + 
 scale_y_continuous(limits = c(-.23, .23)) + 
 scale_x_discrete(labels=c(str_wrap("Venezuela and Colombia share a cultural heritage", width = 20),
              str_wrap("Family history of internal displacement", width = 24),
              str_wrap("Personal history of internal displacement", width = 24),
              str_wrap("Indirect contact with Venezuelans", width = 20),
              str_wrap("Direct contact with Venezuelans", width = 20),
              str_wrap("Cúcuta (vs. Cali)", width = 20),
              str_wrap("Benefits Index", width = 20),
              str_wrap("Salary", width = 20),
              str_wrap("Contract formality", width = 20),
              str_wrap("Skilled labor", width = 20),
              str_wrap("Most Venezuelans in Colombia support guerrillas", width = 25),
              str_wrap("Most Venezuelans in Colombia support the left", width = 25),
              str_wrap("Right Polarization", width = 20)
              )) +
 labs(x = "", 
    y = "Estimate") +
 ggtitle(str_wrap("Openness Index", width = 30)) +
 coord_flip() +
 yy_theme()


### Other questions/indices

## political polarization
# pol_vote1 outcome
pol_vote1_left_ctrls <- lm_robust(update(pol_vote1 ~ pol_left, 
                  reformulate(c(".", labels(demo_ctrls_col)))), data = df_col)

pol_vote1_left_ctrls_plot <- pol_vote1_left_ctrls %>%
    tidy() %>%
    filter(grepl("pol_left", term)) %>%
 ggplot(aes(x = term, y = estimate)) + 
 geom_point(size = 2.5, colour = "#D55E00") + 
 geom_errorbar(aes(ymin = conf.low, ymax = conf.high), 
        width = 0, size = .5, colour = "#D55E00") + 
 geom_hline(aes(yintercept = 0), lty = "dashed") + 
 scale_y_continuous(limits = c(-.5, .5)) + 
 scale_x_discrete(labels=str_wrap("Most Venezuelans in Colombia support the left (1/0)", width = 13)) +
 labs(x = "", 
    y = "Estimate") +
 ggtitle(str_wrap("Most Venezuelans are able to vote in Colombia (1/0)", width = 35)) +
 coord_flip() +
 yy_theme()

## Labor Market index
laborcompet_index_skilled_labor_ctrls <- lm_robust(update(laborcompet_index ~ skilled_labor, 
                  reformulate(c(".", labels(demo_ctrls_col)))), data = df_col)

laborcompet_index_contract_ctrls <- lm_robust(update(laborcompet_index ~ contract, 
                  reformulate(c(".", labels(demo_ctrls_col)))), data = df_col)

laborcompet_index_ctrls_plot <- laborcompet_index_skilled_labor_ctrls %>%
    tidy() %>%
    filter(grepl("skilled_labor", term)) %>%
 bind_rows(
 laborcompet_index_contract_ctrls %>%
    tidy() %>%
    filter(grepl("contract", term))) %>%
 ggplot(aes(x = term, y = estimate)) + 
 geom_point(size = 2.5, colour = "#0072B2") + 
 geom_errorbar(aes(ymin = conf.low, ymax = conf.high), 
        width = 0, size = .5, colour = "#0072B2") + 
 scale_colour_manual(values = c("black")) + 
 geom_hline(aes(yintercept = 0), lty = "dashed") + 
 scale_y_continuous(limits = c(-.5, .5)) + 
 scale_x_discrete(labels=c(str_wrap("Contract formality (1-3)", width = 13),
               str_wrap("Skilled labor (0-2)", width = 13))) +
 labs(x = "", 
    y = "Estimate") +
 ggtitle(str_wrap("Labor Security Index (0-4)", width = 40)) +
 coord_flip() +
 yy_theme()

# Sociotropic
socio_gen_city_ctrls <- lm_robust(update(socio_gen ~ city, 
                  reformulate(c(".", labels(demo_ctrls_col)))), data = df_col)

socio_gen_city_ctrls_plot <- socio_gen_city_ctrls %>%
    tidy() %>%
    filter(grepl("city", term)) %>%
 ggplot(aes(x = term, y = estimate)) + 
 geom_point(size = 2.5, colour = "#000000") + 
 geom_errorbar(aes(ymin = conf.low, ymax = conf.high),
        width = 0, size = .5, colour = "#000000") + 
 scale_colour_manual(values = c("black")) + 
 geom_hline(aes(yintercept = 0), lty = "dashed") + 
 scale_y_continuous(limits = c(-.5, .5)) + 
 scale_x_discrete(labels=str_wrap("Cúcuta (vs. Cali)", width = 13)) +
 labs(x = "", 
    y = "Estimate") +
 ggtitle(str_wrap("Venezuelan migration has been a good thing for the average Colombian (1/0)", width = 40)) +
 coord_flip() +
 yy_theme()

# Cultural index
cultural_index_city_ctrls <- lm_robust(update(cultural_index ~ city, 
                  reformulate(c(".", labels(demo_ctrls_col)))), data = df_col)

cultural_index_city_ctrls_plot <- cultural_index_city_ctrls %>%
    tidy() %>%
    filter(grepl("city", term)) %>%
 ggplot(aes(x = term, y = estimate)) + 
 geom_point(size = 2.5, colour = "#E69F00") + 
 geom_errorbar(aes(ymin = conf.low, ymax = conf.high), 
        width = 0, size = .5, colour = "#E69F00") + 
 scale_colour_manual(values = c("black")) + 
 geom_hline(aes(yintercept = 0), lty = "dashed") + 
 scale_y_continuous(limits = c(-.5, .5)) + 
 scale_x_discrete(labels=str_wrap("Cúcuta (vs. Cali)", width = 13)) +
 labs(x = "", 
    y = "Estimate") +
 ggtitle(str_wrap("Cultural Similarity Index (0-4)", width = 40)) +
 coord_flip() +
 yy_theme()

## plot
open_main_lm.plot 

# +
# (pol_vote1_left_ctrls_plot / laborcompet_index_ctrls_plot /
#  socio_gen_city_ctrls_plot / cultural_index_city_ctrls_plot)

@

Figure \ref{fig:mainreg} presents the correlates of preferences for more open policies toward Venezuelan migrants. Each estimate with 95\% confidence intervals is from a \textit{separate} OLS regression of the openness index regressed on the variable of interest, while controlling for demographic characteristics. For ease of interpretation, we rescale the Openness Index outcome variable and all independent variables shown to vary between 0 and 1.

As is standard in the literature, a respondent's political views are associated with immigration attitudes. Colombians who identify with the political right are less supportive of welcoming border policies, although this does not reach conventional levels of significance. Political misperceptions also are associated with less support for open borders. Colombians who think that Venezuelans support the left or guerrilla groups are less likely to support more welcoming immigration policies. 

We find relatively limited support for common explanations of immigration attitudes. The results do not support the labor market hypothesis, despite the substantial concerns voiced about competition from Venezuelan migrants. A respondent's skill level or labor market contract are not associated with openness towards Venezuelans. We also do not find support for a fiscal burden theory. Respondents' reported salary levels or whether they receive public welfare benefits have little association with their border attitudes. There also is little support for sociotropic effects; respondents in C\'ucuta, which has received a much larger share of Venezuelan migrants, are less likely to say that Venezuelan migration has been good for Colombians.

There is more support for theories based on historical experience, contact, and culture. The coefficient on personal and family histories of displacement are correctly signed, being positively associated with greater openness. Neither reach conventional levels of significance, in part due to our small sample of displaced persons.\footnote{We also excluded Colombian-Venezuelans from the sample as they are citizens of both countries, which reduces those with the most direct experiences.} Colombians who have direct contact with Venezuelans are more likely to support openness than those who only have indirect or no contact. Similarly, those who say that Colombia shares a cultural history with Venezuelan are more likely to support openness. Yet, the presence of large numbers of Venezuelans seems to reify--rather than bridge--differences between groups. Respondents in C\'ucuta, where most Venezuelans have arrived, are less likely to support openness and less likely to think that Venezuelans share a culture than those in Cali (although neither estimate quite reaches conventional levels of significance). Thus, it appears that political fears and direct experiences with Venezuelans through contact are the main correlates with support for tolerant immigration policies. 


%%%%%% REFERENCES %%%%%%%
\newpage
\setstretch{1}
\bibliography{colombia}


\end{document}

